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1.  Introduction 

Satellite-  or  aircraft-borne  synthetic  aperture  radars  (SAR)  have  the  potential  to 
serve  as  a  powerful  and  essential  part  of  the  global  meteorological/oceanographic 
observation  system  (Visecky  and  Stewart,  1982,  Phillips,  1988;  NASA,  1988).  While  the 
potential  of  SAR  systems  is  enormous,  quantitative  interpretation  of  SAR  signals  has 
clearly  been  fiustrated  by  our  incomplete  understanding  of  the  relationships  between  the 
radar  backscatter  cross  section  cf  and  a  complicated  heterogeneous  and  constantly 
changing  state  of  the  sea  surface. 

One  analytical  approach  (Alpers  et  ai,  1981)  is  to  describe  (f  in  terms  of  a 
modulation  transfer  function  /?(K),  in  which  /?  is  a  function  of  wave  tilt  as  well  as 
hydrodynamic  and  orbital  velocity  effects.  It  is  important  to  note,  however,  that  each  of 
these  effects  also  depends  in  a  complicated  fashion  upon  the  structure  of  the  overlying 
marine  atmospheric  boundary  layer  (MABL).  The  use  of  empirical  relationships  such  as 
the  “SASS-1  power  law  model”  has  led  to  considerable  uncertainty  because  the 
fundamental  relationships  are  circumstantially  rather  than  physically  or  mathematically 
defined.  The  combination  of  wave-current  interaction  models  with  physically  based  radar 
scattering  models  is  being  used  by  other  High-Res  Accelerated  Research  Initiative  (ARI) 
Principal  Investigators  to  address  these  problems  Because  the  dominant  atmospheric 
forcing  of  the  sea  surface  is  caused  by  wind  stress,  much  better  knowledge  of  the  spatial 
and  temporal  variability  of  this  stress  is  required  before  the  wave-action  spectrum  can  be 
adequately  specified  by  such  wave-current  interaction  models. 

In  the  first  phase  of  our  High-Res  ARI  work  summarized  here,  we  began 
developing  two  new  marine  atmospheric  boundary  layer  models  of  the  surface  stress 
caused  by  submesoscale  boundary  layer  coherent  structures  and  we  finished  obtaining 
planview  patterns  of  surface  stress  variability  caused  by  MABL  updrafts  and  downdrafts 
(Sikora,  1992,  Sikora  and  Young,  1992,  1993)..  As  a  result  of  observations  taken  in 
September,  1991  during  the  High-Res  ARI  pilot  cruise  in  Gulf  Stream  current,  we  began 
turning  our  attention  to  such  mesoscale  atmospheric  circulations  as  the  solenoidal 
circulation  over  the  sea  surface  temperature  front,  the  coastal  sea  breeze  circulation,  and 
the  flow  between  the  Bermuda  High  and  the  diumally  varying  pressure  trough  on  the 
coastal  plain.  These  phenomena  are  expected  to  occur  again  during  the  second  ARI  cruise 
in  the  summer  of  1993. 

The  submesoscale  convective  structures  of  importance  include  horizontal  roll 
vortices,  which  are  quasi-two-dimensional  circulations  that  fill  the  boundary  layer  (BL), 
and  BL  convective  cells,  which  are  fully  three-dimensional  flows  that  may  replace  the  rolls 
as  the  surface  energy  fluxes  increase,  or  which  may  form  independently  of  the  rolls  in  less 
windy  conditions.  Each  of  these  circulations  depends  sensitively  and  to  varying  degrees 
upon  different  dynamical  and  thermodynamical  forcing  mechanisms;  important  ones 
include  BL  wind  shear,  air/sea  temperature  difference,  and  the  magnitude  of  radiational 
and  latent  heating  in  clouds.  All  of  these  mechanisms  contribute  to  the  transport  of 
horizontal  momentum  downward  and  so  to  stress  variability  at  the  sea  surface  The 
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horizontal  scales  of  this  contribution  to  the  total  stress  are  the  same — about  100  m  to  10 
km — as  that  of  the  circulations  themselves.  Evidence  from  SAR  imagery  of  significant 
stress  variations  caused  by  two-dimensional  BL  roll  vortices  can  be  plausibly  inferred  from 
Fig.  19  of  Visecky  and  Stewart  (1982)  and  is  very  convincingly  found  in  the  pass  1339 
SEASAT  dataset  by  Gerling  (1985,  1986).  Moreover,  three-dimensional  convective  cells 
are  evidenced  in  the  orbit  832  ERS-1  SAR  image  obtained  by  Dr.  Bob  Beal  of  The  Johns 
Hopkins  University  Applied  Research  Laboratory  (APL);  this  image  is  centered  on  the 
coordinates  36°  24'  N,  74°  00'  W  and  was  taken  on  13  September  1991  at  0321  Z  during 
the  first  High-Res  cruise.  Such  submesoscale,  two-  and  three-dimensional  structures, 
which  are  quite  common  over  the  Gulf  Stream,  are  expected  to  be  seen  during  the  1993 
High-Res  field  program  as  well. 

In  this  report,  we  briefly  review  our  progress  on  the  work  that  will  be  continued 
and  extended  during  the  second  phase  of  the  project  from  October  1,  1992  to  September 
30,  1995.  In  Appendix  A  and  Appendix  B  we  give  two  manuscripts  of  journal  articles 
summarizing  our  results.  The  first  one  by  Sikora  and  Young  (1993)  discusses  the 
planview  patterns  of  surface  stress  variability  and  will  appear  in  Boundary  Layer 
Meteorology,  while  the  second  one  by  Wells  et  al.  (1993)  discusses  a  new  method  for 
estimating  the  correlation  dimension  of  boundary  layer  turbulent  time  series. 


2.  Observations  during  the  first  High-Res  cruise 

Instrumentation  for  monitoring  the  surface  fluxes  of  momentum,  heat,  moisture, 
and  downwelling  radiation  was  operated  aboard  the  USNS  Bartlett  during  the  first  High- 
Res  ARI  cruise  in  September,  1991  by  Dr.  George  Young  and  graduate  student  Jeff  Hare. 
The  combination  of  a  sonic  anemometer  for  fast  response  temperature  and  three- 
component  wind  measurements  as  well  as  an  IR  hygrometer  for  fast  response  humidity 
measurements  allowed  computation  of  the  air/sea  fluxes  of  momentum,  heat,  and  moisture 
via  the  inertial  dissipation  method  (Fairall  and  Larsen,  1986).  Jeff  Hare  has  analyzed  these 
observations  under  the  supervision  of  Dr.  Dennis  Thomson,  and  he  reported  his  results  at 
the  spring  1992  Hi-Res  ARI  workshop  that  was  held  at  APL.  These  results  compare  well 
with  measurements  from  a  similar  system  operated  by  Dr.  Jim  Edson  of  Woods  Hole 
Oceanographic  Institute  (WHOI)  aboard  the  RV  Oceanus  during  the  same  experiment. 
However,  owing  to  ship  motions,  observations  of  fluxes  made  on  ships  may  be  hard  to 
interpret,  as  noted  in  the  MS  thesis  by  Hare  (1992). 

Physical  interpretation  of  these  air/sea  flux  variations  across  the  Gulf  Stream  front 
is  complicated  by  mesoscale  atmospheric  evolution  on  the  same  time  scale  as  the  ship’s 
transfrontal  box  patterns.  Thus,  the  meteorological  environment  often  changed 
significantly  in  the  time  required  for  a  ship  to  complete  both  its  frontal  transects  and  its 
along-front  legs  in  the  two  watermasses.  The  resulting  box  of  air/sea  flux  observations 
reflects  a  mixture  of  the  desired  horizontal  gradients  of  air/sea  interaction  with  the 
unavoidable  local  temporal  change  A  primary  accomplishment  during  the  first  phase  of 
this  project  has  been  the  development,  by  Dr.  George  Young,  Todd  Sikora,  and  Jeff  Hare, 
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of  a  method  for  using  the  results  of  the  Penn  State  atmospheric  mesoscale  model,  which 
was  run  in  operational  mode  during  the  experiment,  to  interpret  these  air/sea  flux 
observations.  Further  studies  using  this  mesoscale  model  will  be  performed  in  phase  two 
of  the  project. 

The  signatures  of  several  mesoscale  atmospheric  and  oceanic  phenomena  are 
apparent  in  both  the  observations  and  the  mesoscale  model  results.  The  signatures  in  the 
observations  would  not,  however,  be  interpretable  without  the  vastly  superior  spatial 
coverage  of  the  mesoscale  model.  For  example,  while  the  shipboard  observations  indicate 
that  the  large  amplitude  diurnal  cycle  of  surface  stress  shown  in  Figure  1  is  primarily  a 
consequence  of  a  diurnal  wind  cycle,  they  do  not  indicate  the  origin  of  this  wind  cycle. 
Analysis  of  the  mesoscale  model  results  shows  that  the  wind  cycle  is  a  phase-lagged 
response  to  a  diurnal  cycle  of  the  horizontal  pressure  gradient  in  the  High-Res  region. 
This  diurnal  cycle  in  the  pressure  gradient  between  the  Bermuda  High  and  the  coastal  plain 
is  the  result  of  a  diumally  varying  pressure  trough  over  the  coastal  plain  in  the  lee  of  the 
Appalachian  Mountains.  The  mesoscale  atmospheric  model  consistently  demonstrates  the 
ability  to  correctly  develop  and  simulate  this  temporally  varying  feature  even  when  it  is  not 
adequately  resolved  by  the  synoptic  observations  used  for  model  initialization.  This 
success  stems  from  the  ability  of  the  model  to  predict  the  outcome  of  the  interactions 
between  diumally  varying  solar  heating  and  temporally  unvarying  terrain.  The  model 
typically  requires  about  12  hours  to  overcome  the  errors  in  its  initial  conditions  when  there 
is  such  predictable  diurnal  forcing. 

The  Gulf  Stream  sea  surface  temperature  front  offers  another  mesoscale  forcing 
feature  that  is  essentially  steady  on  the  time  sf'ale  of  the  atmospheric  mesocirculations 
Thus,  it  is  not  surprising  that  the  model  adequately  predicts  the  associated  atmospheric 
mesocirculations  that  were  observed  during  the  first  High-Res  cruise.  Air/sea  flux 
observations  combine  with  the  Penn  State  mesoscale  model  forecasts  and  real-time  clou  J- 
field  observations  from  the  GOES  satellite  and  the  USNS  Bartlett  to  show  that  the  MABL 
responds  vigorously  to  the  sea  surface  temperature  difference  across  the  northwest  all  of 
the  Gulf  Stream.  Significant  differences  in  the  surface  fluxes  across  this  v-atermass 
boundaiy  were  observed  to  result  in  a  sharp  change  in  MABL  structure  across  the 
boundary.  This  atmospheric  front  was  frequently  observed  to  form  in  response  to  the 
oceanic  front  when  benign  synoptic  conditions  brought  atmospheric  flow  along  the 
oceanic  front.  *.ie  resulting  atmospheric  circulation  was  reflected  h\  both  the  enhanced 
cumulus  cloudiness  along  the  front  and  the  convergent  wind  and  stress  patterns  in  the 
mesoscale  model  output.  Analysis  of  this  mesoscale  feedback  between  ocean  and 
atmosphere  is  still  at  a  preliminary  stage.  Use  of  the  mesoscale  model  to  guide 
interpretation  of  the  above  High-Res  observations  will  be  pursued  by  graduate  student 
Sean  Sublette  in  the  second  phase  of  the  project. 
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Figure  1 .  Time  series  of  surface  stress  u.  calculated  using  the  inertial  dissipation  (id)  and 
bulk  aerodynamic  (bulk)  methods  from  observations  taken  on  RV  Oceanus  (O) 
and  Bartlett  (B)  during  the  48-hour  period  beginning  00  EDT  16  September  1991 
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3.  Modeling  of  submesoscale  BL  structures 

In  order  to  identify  the  dominant  forcing  mechanisms  and  characteristic  spatial  and 
temporal  responses  of  the  various  dynamic  and  thermodynamic  modes  of  submesoscale 
MABL  convection,  we  have  conducted  for  a  number  of  years  a  series  of  nonlinear 
modeling  studies  of  BL  rolls  (Shirer,  1986,  Stensrud  and  Shirer,  1988;  Laufersweiler  and 
Shirer,  1989;  Haack  and  Shirer,  1992).  An  important  component  of  these  studies  has  been 
comparison  of  model  results  with  the  MABL  observations  obtained  during  the  KonTur 
and  FIRE  [First  ISSCP  (International  Satellite  Cloud  Climatology  Program)  Regional 
Experiment]  field  programs  (Shirer  and  Brummer,  1986;  Stensrud  and  Shirer,  1988;  Shirer 
and  Haack,  1990).  Emphasis  has  more  recently  been  given  to  studying  the  modification  of 
the  background  wind  and  temperature  profiles  by  the  rolls  in  order  that  the  probable 
forcing  that  produced  the  rolls  could  be  estimated  (Haack  and  Shirer,  1992;  Shirer  and 
Haack,  1990).  These  nonlinear  models  also  yield  profiles  of  the  vertical  fluxes  of 
horizontal  momentum  and  heat.  These  flux  profiles  are  obtained  by  horizontal  integration 
over  one  wavelength  of  the  appropriate  products  given  by  the  model  solutions.  The 
spatial  and  temporal  variability  of  these  products  gives  the  spatial  and  temporal  variability 
of  the  stress  and  vertical  heat  flux  through  application  of  either  standard  bulk  aerodynamic 
or  eddy  correlation  approaches.  These  approaches  are  similar  to  those  used  in  the  analysis 
of  observational  and  Large  Eddy  Simulation  (LES)  results  (e  g.,  Businger  et  al,  1971; 
Moeng,  1984;  Young,  1988) 

The  above  model  results  have  focused  on  the  mid-level  properties  of  the  MABL 
structures,  and  so  it  was  convenient  to  assume  that  the  stress  and  heat  flux  vanished  at  the 
lower  boundary.  These  boundary  conditions  are  clearly  not  suitable  for  study  of  the  stress 
variability  at  the  sea  surface,  and  so  during  the  first  two  years  of  our  ARI  research.  Dr 
Hampton  N.  Shirer,  Dr.  Robert  Wells,  Julie  Schramm,  and  Peter  Bromfield  have  begun 
creating  new  nonlinear,  Boussinesq,  intermediate-order  models  of  both  two-  and  three- 
dimensional  submesoscale  circulations  in  the  MABL.  To  form  the  differential  equations 
for  these  models,  we  first  set  the  bottom  of  the  domain  to  be  at  the  top  of  the  surface 
layer,  assumed  to  be  at  10m  above  the  sea  surface,  and  then  we  choose  lower  boundary 
conditions  that  are  consistent  with  the  wind  and  temperature  profiles  typically  observed  to 
occur  within  the  surface  layer.  Dave  Ledvina  has  written  a  subroutine  that  uses  standard 
similarity  theory  (Liu  et  al.,  1979)  to  give  the  surface  roughness  and  Monin-Obukhov 
length  L  once  the  mean  wind  speed,  air/sea  temperature  difference  and  humidity  difference 
have  been  specified.  We  use  the  values  of  and  L  from  his  subroutine  to  obtain  the 
constants  in  the  boundary  conditions  via  straightforward  manipulations  of  the  similarity 
laws  for  an  unstable  environment  (e  g.,  Stull,  1988).  Other  physical  effects  represented  in 
the  new  models  include  buoyancy,  background  wind  and  temperature  profiles,  Coriolis 
turning,  and  height-dependent  momentum  and  temperature  dissipation. 

To  develop  the  basis  functions  for  the  variable  expansions  used  to  create  the 
intermediate-order  models,  we  solve  eigenfunction  problems  based  on  the  new  boundary 
conditions;  these  problems  yield  nonorthogonal  functions  for  the  vertical  structure  of  the 
variables  that  are  expressible  using  exponential  and  trigonometric  functions  These 
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functions,  as  well  as  the  differential  equations  for  the  models,  are  obtained  symbolically 
with  the  aid  of  the  programs  Derive  and  MAPLE  The  results  are  output  by  Derive  and 
MAPLE  in  FORTRAN  format  for  convenient  insertion  into  the  numerical  codes  for 
simulation  of  the  MABL  circulations.  Sigiuficantly,  this  approach  allows  relatively 
straightforward  extensions  of  the  models.  The  FORTRAN  code  is  in  the  final  stages  of 
debugging;  preliminary  results  indicate  that  the  stress  and  heat  flux  profiles  have  the  same 
qualitative  form  as  that  seen,  for  example,  in  the  submesoscale  MABL  observations  in 
Figure  9  of  Briimmer  and  Busack  (1990) 

However,  because  some  of  the  boundary  conditions  act  as  a  sink  of  energy  while 
others  act  as  a  source,  care  must  be  taken  in  the  way  that  the  perturbation  pressure  is 
handled  in  the  three-dimensional  model.  A  reformulation  of  the  first  version  of  the  model 
created  by  Julie  Schramm  is  an  early  task  of  phase  two  of  the  study.  This  reformulation, 
which  will  be  undertaken  by  graduate  student  Louis  Zuccarello,  will  use  the  vector 
vorticity  equation  as  a  device  for  circumventing  the  need  to  explicitly  represent  the 
pressure  perturbation.  In  addition,  a  separate  two-dimensional  model  based  on  the  along- 
roll  component  of  the  vorticity  is  being  created  by  graduate  student  Peter  Bromfield. 

The  versions  of  the  intermediate-order  models  currently  being  developed  represent 
two-  and  three-dimensional  structures  that  have  a  single  horizontal  wavelength  As  a 
consequence,  the  results  from  these  versions  of  the  model  will  provide  estimates  of  only 
the  symmetric  component  of  sea  surface  stress  variability.  However,  for  the  comparison 
with  observations  using  conditional  sampling  of  updrafts  and  downdrafts  in  the  lower 
portion  of  the  domain,  we  will  need  to  include  at  least  three  horizontal  wavenumbers  in 
the  variable  expansions  for  the  model.  In  this  extended  model  to  be  developed  in  phase 
two  of  the  project,  we  expect  to  maintain  the  same  level  of  resolution  in  the  vertical,  and 
so  this  extension  will  be  relatively  straightforward  to  complete  with  the  aid  of  Derive  or 
MAPLE.  This  new  level  of  resolution  will  be  compatible  with  the  aircraft  data  analysis  by 
Todd  Sikora,  who  gives  the  surface  layer  patterns  of  the  vertical  fluxes  of  heat  and 
momentum  in  updrafts  and  downdrafts  (section  4  and  Appendix  A).  In  our  comparison 
with  these  observations,  we  will  determine  in  particular  whether  the  modeled  stress 
patterns  agree  with  the  observed  ones  that  are  only  weakly  dependent  on  the  degree  of 
instability,  or  thermal  forcing. 

Once  fully  extended,  the  above  models  will  be  ready  for  study  of  both  idealized 
profiles  and  the  ones  observed  during  the  KonTur,  FERE,  ASTEX  and  High-Res  MABL 
field  programs.  To  compare  our  results  with  High-Res  cases,  we  must  supplement  the 
wind  and  temperature  profiles  measured  using  ship-based  sensors  with  profiles  derived 
from  the  Penn  State  mesoscale  model  forecasts,  owing  to  difficulty  in  interpreting  the 
ship-based  measurements  of  these  profiles  (section  2) 

Because  the  initially  developing  solutions  of  the  intermediate-order  model  are 
temporally  periodic  (which  corresponds  to  downwind  propagation  of  the  convective 
circulations),  we  can  obtain  the  expected  wavelength  and  orientation  of  the  MABL 
structures  These  wavelengths  and  orientations  can  be  compared  with  those  inferred  from 
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available  ERS-1  SAR  images  or  visible  satellite  imagery.  The  resulting  solutions  can  be 
conditionally  sampled  in  a  way  similar  to  that  used  by  Todd  Sikora,  who,  during  the  first 
phase  of  the  study,  has  documented  the  coarse  horizontal  variability  of  the  stress  patterns 
in  updrafts  and  downdrafts  measured  at  50  m  above  the  sea  surface  during  FIRE;  these 
results  are  summarized  in  the  next  section  and  Appendix  A.  In  addition,  similar 
conditional  sampling  strategies  will  be  used  to  compare  the  model  solutions  with  the 
sodar-derived  structures  of  the  vertical  velocity  within  thermals  that  were  measured  during 
the  first  High-Res  cruise  and  that  are  expected  to  be  measured  during  the  main  cruise  in 
1993, 


4  Analysis  of  BL  observations 

As  noted  earlier,  mesoscale  circulations  and  submesoscale  MABL  convective 
structures  are  being  emphasized  in  our  observational  analyses  during  phases  one  and  two 
of  our  work  because  they  are  the  primary  atmospheric  modulators  of  surface  stress  on  the 
mesoscale  and  submesoscale  spatial  and  temporal  scales  in  the  High-Res  region  of  the 
Gulf  Stream  The  primary  goal  of  our  analyses  is  the  development  of  quantitative 
conceptual  models  of  these  mesoscale  and  submesoscale  circulations  within  different 
large-scale  atmospheric  regimes.  The  major  objective  is  the  determination  of  the  temporal 
and  spatial  patterns  of  mesoscale  and  submesoscale  surface  stress  and  heat  flux  variability 
for  each  observed  case,  stratified  where  possible  by  appropriate  large-scale  forcing  as 
measured  by  the  controlling  parameters. 

The  quantitative  observational  description  of  the  role  of  submesoscale  MABL 
convection  in  causing  horizontal  air/sea  flux  variability  was  summarized  in  the  MS  thesis 
by  Sikora  (1992),  was  presented  at  the  Tenth  Symposium  on  Turbulence  and  Diffusion 
(Sikora  and  Young,  1992)  and  will  appear  in  Boundary  Layer  Meteorology  (Sikora  and 
Young,  1993,  which  is  Appendix  A)  Aircraft  turbulence  measurements  made  during  the 
first  stratocumulus  phase  of  FIRE  (Albrecht  et  al ,  1988)  that  was  conducted  under  the 
1986-1991  ONR-fiinded  University  Research  Initiative  (URI)  at  Penn  State  provided  the 
primary  data  source  Todd  Sikora  in  his  MS  research  (Sikora,  1992)  prepared  horizontal 
planviews  depicting  the  typical  surface  layer  patterns  of  the  vertical  fluxes  of  momentum, 
heat,  and  moisture  associated  with  the  updraft  and  downdraft  components  of 
submesoscale  convective  structures  in  the  MABL.  For  example.  Figure  2  shows  the 
corresponding  composite  planview  perturbation  wind  field  for  MABL  convective 
downdrafts,  the  phenomenon  that  causes  "cat’s  paw"  wave  patterns  at  the  sea  surface 
Physical  interpretation  of  these  flux  patterns  shows  that  these  flux  patterns  exhibit  only 
limited  sensitivity  to  large-scale  atmospheric  and  oceanic  conditions.  While  the  patterns 
can  differ  significantly  between  stable  and  unstable  atmospheric  boundary  layers,  they  are 
essentially  similar  over  a  broad  range  of  unstable  conditions  Thus,  these  submesoscale 
flux  patterns  are  expected  to  be  rather  ubiquitous  within  the  MABL  These  results  can  be 
used  by  ocean-wave  modelers  to  quantify  the  horizontal  variability  of  atmospheric  forcing 
on  spatial  scales  on  the  order  of  0  1  to  1  0  km. 
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Plaoview  ©f  Composite 
Dowodiraft  Horizootai  Wind 


Figure  2.  Planview  of  composite  horizontal  wind  the  MABL  convective  downdrafts  in 
the  surface  layer  (after  Sikora  and  Young,  1993). 
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As  the  observational  study  of  MABL  convection  is  complete,  the  only  task 
remaining  in  this  aspect  of  the  project  is  to  use  these  results  in  the  modeling  effort 
described  above  in  section  3.  During  the  second  phase  of  this  work,  these  results  will  be 
used  to  guide  the  final  stages  of  model  development  as  well  as  for  the  verification  and 
interpretation  of  model  results.  Of  particular  interest  is  whether  the  model  corroborates 
the  observed  relatively  weak  dependence  of  the  flux  patterns  on  the  degree  of  instability. 


5.  An  improved  correlation  dimension  algorithm 

The  atmospheric  structures  responsible  for  creating  stress  variability  at  the  sea 
surface  are  clearly  chaotic,  and  there  remains  a  major  challenge  for  scientists  to  suitably 
quantify  this  chaotic  structure.  Existing  algorithms  for  estimating  the  correlation 
dimension  d  typically  produce  results  that  are  difficult  to  interpret  and  that  are  subject  to 
errors  whose  magnitudes  are  not  typically  estimated.  We  have  developed  a  new  objective 
algorithm  (Wells  et  al.,  1993  and  Appendix  B)  that  overcomes  these  difficulties  because  it 
is  based  on  hypotheses  that  are  more  likely  to  be  satisfied  than  are  the  ones  for  the 
standard  approach.  Not  only  does  this  method  in  principle  produce  an  infinite  number  of 
estimates  for  the  correlation  dimension,  but  it  also  gives  error  estimates  that,  when 
minimized,  produce  the  optimal  value  for  d.  This  algorithm  is  developed  and  tested  in  a 
manuscript  that  has  been  submitted  to  Physica  D  (Wells  et  al.,  1993  and  Appendix  B),  a 
preliminary  version  was  presented  as  a  poster  at  the  Spring  1992  AGU  meeting  in 
Montreal  (Wells  et  a/.,  1992).  Christian  Fosmire  in  his  MS  thesis  work  (Fosmire,  1993) 
tests  this  algorithm  using  time  series  data  from  both  the  standard  Lorenz  and  Henon 
attractors  and  from  sodar  data  measured  in  a  1988  field  experiment  conducted  near  Penn 
State.  In  addition,  in  collaboration  with  Dr  Dennis  Thomson,  Dr  Harry  Henderson  has 
been  applying  this  and  other  algorithms  to  BL  measurements  (eg.,  Thomson  and 
Henderson,  1991,  1992;  Thomson  ^ro/.,  1991). 
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1.  Introduction 

The  marine  atmospheric  surface  layer  (MASL)  (Z  <  O.l 
where  is  the  boundary  layer  depth)  (Stull,  1989)  plays  a 
major  role  in  air/sea  interaction.  It  is  the  MASL  that 
couples  the  sea  surface  to  the  overlying  atmospheric  boundary 
layer  through  vertical  eddy  fluxes  (hereafter  referred  to 
simply  as  fluxes) ,  such  as  those  of  buoyancy,  moisture,  and 
momentum.  While  these  fluxes  are  present  in  both  stable  and 
unstable  environments,  they  tend  to  be  enhanced  by  convective 
elements,  as  is  implied  by  bulk  aerodynamic  parameterization 
(e.g.  Liu  et  al.,  1979).  For  this  reason,  much  of  the  recent 
research  in  this  area  has  concentrated  on  the  convective 
marine  atmospheric  surface  layer  (CMASL) . 

CMASL  fluxes  are  realized  through  sub-mesoscale 
convective  updrafts  (CUs)  and  downdrafts  (CDs) .  These 
features  have  diameters  on  the  order  of  tens  to  thousands  of 
meters  (Lenschow  and  Stephens,  1980) .  Khalsa  and  Greenhut 
(1985)  found  that  within  a  central  Pacific  CMASL,  such 
features  were  responsible  for  75%  of  the  total  flux  of  heat, 
moisture,  and  momentum.  This  imp.rtance  provides  the 
motivation  to  explore  further  the  flux  characteristics  of 
CMASL  drafts. 

Various  methods  may  be  employed  to  study  CMASL  drafts. 
Qualitatively,  this  can  be  accomplished  by  simply  observing 
the  effects  their  associated  fluxes  have  on  the  environment. 
On  a  human  scale,  anyone  who  has  enjoyed  a  day  of  sailing, 
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seen  wisps  of  sea  fog,  or  witnessed  cats  paws  rippling  across 
a  body  of  water  has  personally  sensed  the  impact  of  CMASL 
drafts  on  the  environment.  Another  example  of  this 
manifestation  can  be  seen  in  Figures  1  and  2  of  Sikora  and 
Young  (submitted  to  Boundary  Layer  Meteorology,  November 
1992)  .  Figure  1  is  a  visible  image  from  the  Kosmos  1500 
(Okean)  satellite.  The  image  is  centered  at  26°  N  125.5°  W 
and  is  dated  July  11,  1984.  It  shows  common  boundary  layer 
cumuliform  clouds,  horizontal  scales  of  which  are  on  order  of 
the  boundary  layer  depth,  induced  by  fluxes  of  heat  and 
moisture  within  CMASL-born  CUs.  Figure  2  is  the  corresponding 
X-band  real  aperture  radar  image  of  the  underlying  sea 
surface.  It  depicts  a  perturbed  sea  surface  wave  pattern 
resulting  from  CMASL  momentum  flux  patterns  driven  by  flow 
into  CUs  and  out  of  CDs. 

In  order  to  obtain  data  on  the  effects  and  structure  of 
CUs  and  CDs  in  a  more  quantitative  sense,  various  techniques 
have  been  employed,  including  observations  from  a  sensor 
equipped  catamaran  (Dorman  and  Mollo-Christensen,  1973) ,  the 
use  of  acoustic  sounders  (Gaynor  and  Mandics,  1978) ,  and 
observations  of  sea  gulls  in  flight  (Woodcock,  1940,  1975) . 
Data  for  the  quantitative  study  in  this  paper  are  obtained 
from  NCAR  Electra  aircraft  flights  during  the  First  ISCCP 
(International  Satellite  Cloud  Climatology  Project)  Regional 
Experiment  (FIRE) .  Flights  were  conducted  over  a  several 
hundred  kilometer  region  west  of  the  southern  California  coast 
during  June  and  July,  1987.  For  detailed  reviews  of  the  FIRE 
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project,  see  Albrecht  et  al.  (1988)  and  Kloesel  et  al.  (1988)  . 

The  traditional  way  to  investigate  the  flux 
characteristics  of  boundary  layer  features  is  through 
conditional  sampling  and  composite  analysis,  examples  of  which 
are  discussed  below.  Simply  put,  conditional  sampling  uses  an 
indicator  function  to  isolate  from  a  dataset  features  of 
interest  (e.g.  CUs  and  CDs) .  Composite  analysis  then  combines 
information  from  these  selected  features  to  determine  their 
average  spatial  structure. 

In  the  past,  various  studies  have  gathered  valuable 
information  concerning  boundary  layer  features  using  these 
sampling  and  analysis  technigues.  Wilczak  and  Tillman  (1980) 
and  Wilczak  (1984)  used  Boulder  Atmospheric  Observatory  tower 
data  to  produce  detailed,  land-based,  statistics  for 
convective  atmospheric  surface  layer  features  known  as 
temperature  ramps.  Young  (1988)  also  investigated  land-based 
turbulence  at  the  Boulder  Atmospheric  Observatory.  He 
conditionally  sampled  aircraft  data,  based  upon  vertical 
velocity  perturbations  and  mixed  layer  spectra  of  vertical 
velocity  and  temperature,  in  order  to  analyze  thermals  within 
the  convective  boundary  layer.  Lenschow  and  Stephens  (1980) 
used  NCAK  Electra  aircraft  data  to  generate  information  on  the 
structure  of  marine  atmospheric  boundary  layer  (MABL)  thermals 
using  humidity  as  an  indicator  function.  Greenhut  and  Khalsa 
(1982)  and  Khalsa  and  Greenhut  (1985)  also  used  aircraft  data 
to  study  properties  of  drafts  in  the  MABL  but,  as  in  Young 
(1988),  they  used  vertical  velocity  as  an  indicator  of  a 
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feature.  In  addition  to  observational  studies,  large-eddy 
simulation  models,  such  as  that  of  Schumann  and  Moeng  (1991), 
have  been  used  to  provide  boundary  layer  datasets  for 
conditional  sampling  and  compositing. 

While  previous  work  has  resulted  in  detailed  vertical 
profiles  of  properties  of  MABL  updrafts  and  downdrafts, 
information  on  planview  flux  variability  of  CMASL  CUs  and  CDs 
is  lacking.  To  fill  this  need,  this  study  utilizes 
conditionally  sampled  and  composited  eddy  correlated  FIRE 
aircraft  data  to  investigate  the  typical  CMASL  CU  and  CD 
planview  patterns  of  vertical  velocity  flux,  w'w' ,  buoyancy 
flux,  w'T^',  absolute  humidity  flux,  w'r',  and  along-mean-wind 
momentum  flux,  w'u'.  Note  that  while  CUs  and  CDs  are 
segregated  in  this  study,  the  technigues  used  for  their 
analyses  are  identical. 
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2 .  Procedures 

The  use  of  NCAR  Electra  turbulence  data  from  project  FIRE 
is  advantageous  for  the  type  of  research  discussed  here 
because  a  properly  instrumented  aircraft  provides  a  means  for 
collecting  high  resolution  turbulence  data  from  a  large  number 
of  CMASL  drafts.  While  each  draft  can  be  penetrated  only 
once,  penetrations  of  many  different  drafts  can  be  made  from 
numerous  angles.  In  addition,  the  turbulence  data  from  a 
draft  sampled  by  an  aircraft  can  be  looked  upon  as  a  snapshot 
of  that  feature  (see  sub-section  2.3.).  The  result,  after 
data  processing,  conditional  sampling,  and  compositing,  is  a 
statistically  robust  planview  flux  description  of  a  typical 
CMASL  draft. 

2.1.  FLIGHT  INFORMATION 

Of  the  numerous  available  datasets  collected  from  Electra 
flight  legs  during  FIRE,  only  the  20  used  in  this  study  were 
conducted  within  the  CMASL.  These  legs  averaged  50  km  (500 
sec)  in  length  (duration)  .  These  legs  were  flown  on  seven 
different  days,  during  which  the  average  boundary  layer  depth 
was  920  m.  The  legs  were  flown  at  a  height  of  50  m,  well 
within  the  depth  of  the  surface  layer.  The  ratio  of  the 
boundary  layer  depths  to  the  Monin-Obukhov  lengths  ranged  from 
-0.62  to  -37.56,  indicative  of  slightly  to  moderately  unstable 
MABLs.  While  flux  intensity  differences  for  CUs  and  CDs  are 
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seen  within  this  stability  range,  the  "flux-shapes"  (the 
patterns  of  spatial  variation  of  flux  magnitude  across  a 
draft)  are  similar.  This  pattern  similarity  in  the  presence 
of  intensity  differences  is  exactly  analogous  the  geometric 
similarity  of  "similar  triangles"  of  varying  sizes. 
Structural  differences  are,  however,  observed  for  stable  and 
extremely  unstable  cases.  Insufficient  data  prevented  the 
extension  of  the  composite  analysis  to  these  other  stability 
ranges . 

2.2.  INSTRUMENTATION  AND  DATA  PROCESSING 

The  Electra  CMASL  20  hz  turbulence  data  used  in  the 
present  study  include  vertical  velocity,  w,  air  temperature, 
T,  absolute  humidity,  r,  and  the  aircraft-oriented  components 
of  the  wind,  u,  and,  v.  This  sampling  rate  combined  with  the 
100  m/s  aircraft  speed  yields  a  minimum  resolvable  wavelength 
of  10  m.  A  review  of  the  aircraft  instrumentation  of 
significance  to  this  study  can  be  found  in  Nucciarone  and 
Young  (19S1)  .  Data  processing  for  mean,  trend,  and  spike 
removal,  as  well  the  derivation  of  buoyancy  flux,  parallels 
that  of  Moyer  and  Young  (1991)  .  In  order  to  eliminate 
contributions  Ly  mesoscale  phenomena  in  the  data,  high-pass 
filtering  is  employed  following  the  techniques  of  Young 
(1987,1988) . 

An  example  of  the  spectral  response  of  the  filter  can  be 
seen  in  Figure  1.  The  figure  shows,  for  a  cosine  wave,  how 
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Fig.  1.  An  example  of  the  spectral  response  of  the  filter. 

the  ratio  of  the  wave's  filtered  variance  to  its  unfiltered 
variance  varies  with  the  ratio  of  wavelength  to  In  this 
example,  the  cutoff  wavelength  chosen  is  0.1  Z^.  At  this 
wavelength,  50%  of  the  wave's  original  variance  remains.  At 
longer  wavelengths,  the  amount  of  the  wave's  original  variance 
remaining  after  filtering  diminishes  rapidly. 

FIRE  spectra  of  Nucciarone  and  Young  (1991)  indicate  that 
1.5  is  the  appropriate  mesoscale  subrange  cutoff  wavelength 
for  use  in  the  filter.  The  net  result  of  this  data 
processing,  for  any  one  variable,  is  a  sub-mesoscale 
perturbation  series  (denoted  by  primed  variables)  suitable  for 
computing  eddy  correlation  statistics. 

Coordinate  rotations  are  performed  on  the  perturbation 
wind  components  yielding  one  component  in  the  direction  of  the 
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mean  wind,  ,  and  another  aligned  with  positive  values  90°  to 
the  left  of  the  mean  wind,  v' .  Unlike  that  seen  in  Geernaert 
and  Hansan  (1992),  where  orographic  distortion  of  the  low 
level  flow  may  be  occurring,  the  cross-mean-wind  components  of 
the  momentum  fluxes  of  this  study  are  both  statistically  and 
physically  insignificant  when  compared  with  the  along-mean- 
wind  component  (i.e.  mean  <  standard  deviation  for  w'v'  and 
w'v'  «  w'u').  For  this  reason,  the  cross-mean-wind  momentum 
flux  will  not  be  discussed. 

2.3.  CONDITIONAL  SAMPLING 

In  order  to  distinguish  updrafts  and  downdrafts,  event 
criteria  based  on  perturbation  vertical  velocity  are  required. 
Greenhut  and  Khalsa  (1982)  and  Young  (1988)  show  that  w'  event 
criteria,  incorporating  germane  horizontal  scales  and 
magnitude  thresholds,  allow  for  the  proper  detection  of 
boundary  layer  features  of  interest. 

In  order  to  distinguish  CUs  and  CDs,  the  w'  series  is 
first  band-pass  filtered  to  retain  only  the  convective  scales 
(Young  1988) .  The  band-pass  filter  is  a  variation  on  that 
discussed  in  sub-section  2.2.,  combining  both  a  high-pass  and 
low-pass  stage.  It  is  designed  to  eliminate  both  the 
mesoscale  and  inertial  subrange  contributions  to  the  w' 
series,  while  still  preserving  the  CUs  and  CDs  of  the  energy 
containing  subrange.  As  discussed  in  sub-section  2.2.,  1.5 
is  used  as  the  mesoscale  subrange  cutoff  wavelength.  The 
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choice  of  an  inertial  subrange  cutoff  wavelength  is  0.1  Z^. 
As  in  Young  (1988) ,  this  choice  eliminates  scales  of  motion 
associated  with  the  inertial  cascade. 

Given  the  average  of  this  study,  a  minimum  draft  width 
of  50  m  is  chosen  so  that  the  narrowest  acceptable  event 
(CU/CD  couplet)  corresponds  with  the  0.1  inertial  subrange 
cutoff  wavelength  used  in  the  band-pass  filter.  All  drafts 
not  meeting  the  minimum  width  criteria  are  rejected  because  of 
their  inertial  subrange  character. 

The  literature  mentioned  in  section  1.  provides  numerous 
examples  of  conditional  sampling  magnitude  thresholds  designed 
to  distinguish  features  of  interest  from  their  environment. 
Inspection  of  the  FIRE  w'  data  series  suggests  that  a  w' 
threshold  of  ±  0.1  m/s  eliminates  disorganized  areas  of  weak 
ascent  or  decent.  Any  data  points  for  which  w'  <  0.1  m/s  are 
therefore  considered  not  to  be  part  of  significant  convective 
features  and  are  so  rejected. 

CUs  (CDs)  are  then  defined  to  be  features  within  the 
band-pass  filtered  w'  data  series  meeting  these  minimum  width 
and  magnitude  criteria.  The  total  number  of  CUs  sampled  and 
usev.>  in  compositing  is  2839,  averaging  98  m  in  width  and 
occupying  28%  of  the  data  series.  The  total  number  of  CDs 
sampled  and  used  in  compositing  is  3062,  averaging  107  m  in 
width  and  occupying  33%  of  the  data  series.  Given  the  above 
widths  and  an  aircraft  velocity  of  100  m/s,  it  can  be  shown 
that  Taylor's  Hypothesis  is  valid  (Stull  1989)  .  Draft 
composites  of  eddy  correlation  statistics  can  therefore  be 
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looked  upon  as  typical  planview  flux  snapshots  of  a  CU  or  a 
CD. 


2.4.  COMPOSITING  OF  EDDY  CORRELATION  STATISTICS 

After  the  isolation  of  CUs  and  CDs,  eddy  correlation 
fluxes  of  vertical  velocity,  buoyancy,  absolute  humidity,  and 
along-mean-wind  momentum  are  calculated  from  the  high-pass 
filtered  dataset,  for  each  draft.  Note  that,  unlike  the  band- 
passed  series  used  to  locate  the  sub-mesoscale  CUs  and  CDs, 
this  series  retains  the  inertial  subrange  contributions.  Only 
the  mesoscale  contributions  are  filtered  out,  as  discussed  in 
sub-section  2.2..  In  order  to  develop  typical  planview  flux 
snapshots,  the  fluxes  from  all  the  drafts  are  composited 
together  using  the  series  of  averaging  calculations  and 
grouping  techniques  discussed  below. 

Each  draft  is  first  divided  into  3  bins  of  equal  length: 
the  aircraft  entry  region  of  the  draft,  bin  (a),  the  middle  of 
the  draft,  bin  (b)  ,  and  the  aircraft  exit  region  of  the  draft, 
bin  (c)  .  The  average  CU  and  CD  flux  patterns  for  each  leg  are 
then  found  ^.y  calculating  bin  average  eddy  correlation 
statistics  as  follows.  The  flux  of  any  variable  y  by  variable 
w'  is  found  by  first  taking  the  sum  of  their  product  over  all 
points  of  all  like  bins  [bin  (a),  (b) ,  or  (c) ]  over  all  like 
drafts  (CUs  or  CDs)  within  the  leg.  This  sum  is  then  divided 
by  the  number  of  data  points  in  all  like  bins  of  all  like 
drafts  within  the  leg. 
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The  resulting  40  leg  average  drafts  (LADs)  from  the  20 
flight  legs  are  then  separated  into  one  of  3  equal  angular- 
width  groups.  The  grouping  of  an  LAD  depends  on  the  path 
taken  by  the  aircraft,  through  the  draft,  relative  to  the  mean 
wind.  In  polar  coordinates  with  the  mean  wind  vector  aligned 
towards  0°,  these  relative  heading  groups  are:  along-mean- 
wind  (330°  to  30°  and  150°  to  210°),  cross-mean-wind  (60°  to 
120°  and  240°  to  300°) ,  and  diagonal-to-the-mean-wind  (the 
remaining  relative  heading  ranges) . 

To  maintain  wind-relative  consistency  in  bin  naming,  flux 
data  for  bins  (a)  and  (c)  are  switched  for  those  drafts  of  the 
along-mean-wind  and  diagonal  groups  whose  relative  headings 
oppose  the  mean  wind.  The  same  is  done  for  LADs  of  the  cross- 
mean-wind  group  whose  relative  heading  range  is  aligned  to  the 
left  of  the  mean  wind.  Note  that  while  there  are  twice  as 
many  relative  heading  ranges  for  the  diagonal  group  as  for  the 
other  two,  the  angular  width  of  each  heading  range  in  the 
diagonal  group  is  only  half  as  large  as  for  the  other  two. 
Thus,  the  angular  coverage  of  all  3  groups  is  equal. 

The  two  remaining  relative-heading  ranges  of  the  diagonal 
group  (30°  to  60°  aiid  300°  to  330°)  can  be  merged  using 
symmetry  arguments.  Wilczak  (1984)  shows  mirror  symmetry  of 
the  horizontal  perturbation  wind  field  about  an  axis  aligned 
with  the  mean  wind  extending  through  the  center  of  his  large- 
scale  eddy  (LSE  which  is  a  CU/CD  couplet) .  A  similar  symmetry 
can  be  seen  in  the  current  dataset  by  examination  of  the 
fluxes  of  the  cross-mean-wind  group  composites  found  in 
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section  3 . .  Comparison  of  the  data  from  the  two  relative 
heading  ranges  of  the  diagonal  group  (not  shown)  demonstrates 
symmetry  across  the  axis  of  the  mean  wind  for  both  the  upwind 
and  downwind  ends  of  the  drafts.  Taking  these  symmetries  into 
account,  flux  data  of  like  LAD  downwind  bins  of  the  diagonal 
group  are  composited  together  as  are  flux  data  of  upwind  bins. 
Figure  2  shows  the  spatial  relation  of  the  three  groups 
within  a  planview  draft.  Bins  are  labeled  (a) ,  (b) ,  or  (c) . 
Upper  case  letters  indicate  the  group  to  which  each  bin 
belongs;  the  a long-mean-wind  group  (A)  ,  the  cross-mean-wind 
group  (C) ,  and  the  diagonal-to-the-mean-wind  group  (D) .  The 
mean  wind  vector  is  assumed  to  be  directed  towards  the  top  of 
the  figure. 

The  process  of  compositing  flux  data  of  like  LAD  bins 
into  a  flux  for  the  corresponding  bins  of  a  group  average 
draft  (GAD)  will  now  be  explained.  First,  because  the  center 
bin  for  all  groups  corresponds  to  the  same  part  of  the  draft, 
the  average  bin  (b)  flux  over  all  like  LADs  over  all  groups 
(bin  (b)  average)  may  be  used  to  represent  the  flux  in  the 
center  bin  of  all  groups.  As  stated  in  sub-section  2.1., 
drafts  of  each  leg  are  structurdly  similar  in  terms  of  their 
flux-shape.  However,  there  is  a  need  to  correct  for  mean 
draft  intensity  differences  between  like  bins  of  different 
LADs.  It  follows  then  that  LAD  bin  (a)  and  (c)  fluxes  are 
rescaled  before  any  averaging  to  produce  GADs.  This  rescaling 
uses  a  correction  coefficient  based  on  the  need  for  inter-leg 
similarity  in  bin  (b)  ,  described  above.  This  correction 
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coefficient  is  bin  (b)  average  divided  by  the  respective  LAD 
bin  (b)  .  This  approach  is  robust  because  the  LAD  bin  (b) 
fluxes  of  all  the  quantities  studied  are  generally  large  and 
of  same  sign. 

The  final  planview  composite  is  found  by  simply 
overlaying  the  along-mean-wind,  cross-mean-wind,  and  diagonal- 
to-the-mean  wind  (again  using  symmetry  for  the  second 
diagonal)  GADs  for  each  statistic  discussed. 
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Fig.  2.  The  spatial  relation  of  the  three  groups  within  a 

planview  draft. 
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3 .  Observational  Results 

Planview  draft  flux  composites  of  vertical  velocity, 
buoyancy,  absolute  humidity,  and  the  a  long-mean-wind  component 
of  momentum,  along  with  corresponding  bin  standard  deviations, 
are  presented  in  Tables  I  through  VIII.  Each  table  is  a  three 
by  three  array  of  group  bin  fluxes.  Bin  location  in  this 
array  corresponds  to  its  position  within  the  composite  draft. 
In  all  tables,  the  mean  wind  is  directed  towards  the  top  of 
the  table.  For  each  composite,  then,  the  along-mean-wind 
group  axis  is  aligned  from  the  upwind  center  to  downwind 
center,  the  cross-mean-wind  group  axis  is  aligned  from  the 
middle  left  to  middle  right,  and  the  diagonal  group  makes  up 
the  remainder. 

Note  that  the  middle  center  bin  of  all  tables  is  its 
respective  bin  (b)  average  of  sub-section  2.4..  Recall  also 
that  the  middle  center  bin  is  used  to  rescale  all  other  bins 
within  a  flux  composite.  For  these  reasons,  the  majority  of 
the  uncertainty  associated  with  any  flux-shape  and  intensity 
related  inter leg  differences,  associated  with  inter leg 
differences  in  the  middle  center  bins,  is  forctU  into  the 
perimeter  of  the  composite  array.  It  follows  then  that 
standard  deviations  for  the  middle  center  bins  are  not  given. 


3.1.  VERTICAL  VELOCITY  FLUX 


Vertical  velocity  dictates,  to  a  large  extent,  the  other 
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flux  patterns  discussed  in  this  study.  Not  only  is  it  used 
for  the  defining  of  events,  but  also  it  is  the  flux  variable 
doing  the  work,  so  to  speak,  by  advecting  the  other 
quantities.  The  planview  CU  w'w'  composite  is  presented  in 
Table  I. 


TABLE  I 

CU  vertical  velocity  flux  composite  and 
corresponding  standard  deviations. 


w'w' (m^s 
(std.  dev.) 

Left 

Center 

Right 

Downwind 

0.227 

(0.015) 

0.244 

(0.040) 

0.227 

(0.015) 

Middle 

0.217 

(0.023) 

0.341 

0.225 

(0.021) 

Upwind 

0.238 

(0.009) 

0.266 

(0.016) 

0.238 

(0.009) 

Because  all  bin  standard  deviations  are  much  smaller  than 
the  bin  fluxes,  there  is  statistical  confidence  in  the  entire 
CU  w'w'  composite.  The  strongest  flux  region  is  found  along 
the  along-mean-wind  group  axis  with  the  largest  magnitude  of 
flux  being  at  the  center  of  the  composite.  This  shape 
compares  well  with  that  found  at  50  m  over  land  by  Wilczak 
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(1984) .  Examining  the  perimeter  of  the  composite,  the  flux 
magnitude  along  the  upwind  side  of  the  draft  tends  to  be 
larger  than  that  along  the  downwind  side,  although  these 
asymmetries  are  less  statistically  significant  than  the  radial 
variation.  The  weakest  flux  is  found  in  the  middle  right  and 
left  bins. 

The  planview  CD  vertical  velocity  flux  composite  is 
shown  in  Table  II.  As  with  the  CU  composite,  there  is 


TABLE  II 

CD  vertical  velocity  flux  composite  and 
corresponding  standard  deviations. 


w'w'  (m^s”^) 

( std .  dev . ) 

Left 

Center 

Right 

Downwind 

0.160 

(0.011) 

0.161 

(0.018) 

0.160 

(0.011) 

Middle 

0.171 

(0.024) 

0.204 

0.173 

(0.024) 

Upwind 

0.162 

(0.009) 

0.166 

(0.019) 

0.162 

(0.009) 

statistical  confidence  in  all  bin  fluxes.  While  the  largest 
magnitude  of  w'w'  is  again  found  in  the  center  middle  bin,  the 
over  all  flux  composite  is  weaker  than  that  of  the  CU.  This 
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finding  is  in  agreement  with  those  of  Greenhut  and  Khalsa 
(1982)  for  their  coherent  downdrafts.  As  with  the  CU 
composite,  other  aspects  of  the  w'w'  pattern  are  less 
statistically  significant.  The  strongest  flux  region  is  found 
along  the  cross-mean-wind  group  axis.  Flux  along  the  upwind 
side  of  the  draft  is  slightly  larger  than  that  along  the 
downwind  side. 

3.2.  BUOYANCY  FLUX 


The  buoyancy  flux  of  the  current  study  has  the  same  sign 
and  approximates  to  a  large  extent  the  magnitude  of  the  heat- 
flux  w'T' .  The  planview  cu  buoyancy  flux  composite  is 
presented  in  Table  III.  While  the  composite  corners  lack 
statistical  significance,  there  is  statistical  confidence  in 
that  part  of  the  composite  consisting  of  the  along-mean-wind 
and  cross-mean-wind  groups.  Along  these  axes,  general 
symmetries  exist  with  strongest  w'T^'  being  at  the  middle 
center  of  the  composite.  All  fluxes  within  the  composite  are 
down-gradient.  This  finding  is  in  agreement  with  those  of 
Khalsa  and  Greenhut  (1985)  for  their  lowest  levels. 

The  planview  CD  buoyancy  flux  composite  is  presented  in 
Table  IV.  Statistical  significance  in  the  w'T^'  pattern  is 
similar  to  that  of  the  cu.  The  along-mean-wind  and  cross¬ 
mean-wind  group  axes  are  more  symmetric  than  in  the  CU  while 
the  diagonal-to-the-mean-wind  group  axis  is  even  more 
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TABLE  III 

CU  buoyancy  flux  composite  and  corresponding 
standard  deviations. 


w'T^'  (ms"^°C) 
(std.  dev.) 

Left 

Center 

Right 

Downwind 

0.024 

(0.021) 

0.010 

(0.002) 

0.024 

(0.021) 

Middle 

0.010 

(0.002) 

0.018 

0.011 

(0.002) 

Upwind 

0.006 

(0.008) 

0.014 

(0.002) 

0.006 

(0.008) 

asymmetric  in  the  along-mean-wind  direction  than  in  the  CU. 

The  pattern  of  the  CD  w'T^'  magnitude  is  similar  to  that 
of  the  CU.  However,  not  all  fluxes  are  positive.  Fluxes 
found  within  the  upwind  left  and  right  bins,  although  not 
statistically  significant,  are  counter-gradient.  In  contrast 
and  as  will  be  seen  in  sub-section  3.3.,  all  fluxes  composing 
the  CD  absolute  humidity  composite  CU  of  this  study  are  down- 
gradient.  Keeping  in  mind  the  lack  of  statistical  confidence, 
the  upwind  left  and  right  bins  of  w'T^'  for  the  CD  composite, 
then,  depart  from  findings  of  Khalsa  and  Greenhut  (1985)  that 
warm  dry  downdrafts  have  negative  buoyancy  flux  only  at  levels 
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TABLE  IV 

CD  buoyancy  flux  composite  and  corresponding 
standard  deviations. 


w'Ty'  (ms"^°C) 

( std .  dev . ) 

Left 

Center 

Right 

Downwind 

0.009 

(0.005) 

0.005 

(0.001) 

0.009 

(0.005) 

Middle 

0.004 

(0.001) 

0.006 

0.004 

(0.001) 

Upwind 

-0.002 

(0.009) 

0.004 

(0.001) 

-0.002 

(0.009) 

higher  than  the  surface  layer.  Similar  to  what  is  suggested 
in  Khalsa  and  Greenhut  (1985)  for  some  of  their  anomalous 
patterns,  the  counter-gradient  of  this  study  may  result 
from  entrainment  of  buoyant  CU  air  into  the  upwind  side  of  the 
CD.  The  positive  buoyancy  flux  portion  of  the  CD  w'T^' 
composite  of  this  study  compares  well  with  findings  in  the 
lowest  levels  of  Khalsa  and  Greenhut  (1985) . 

Of  final  note:  a  surface  layer  counter-gradient  heat 
flux  is  found  in  the  upwind  region  of  the  ensemble  downdraft 
of  Wilczak's  (1984)  vertical  cross  section.  Thus,  the 
finding  of  Wilczak  (1984)  compares  well  with  the  current 
study . 
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3.3.  ABSOLUTE  HUMIDITY  FLUX 


The  planview  CU  absolute  humidity  flux  composite  is 
presented  in  Table  V.  Statistical  confidence  in  w'r'  exists 


TABLE  V 

CU  absolute  humidity  flux  composite  and 
corresponding  standard  deviations. 


w'r'  (gs~^m“^) 

( std .  dev . ) 

Left 

Center 

Right 

Downwind 

0.045 

(0.003) 

0.041 

(0.005) 

0.045 

(0.003) 

Middle 

0.039 

(0.005) 

0.068 

0.041 

(0.008) 

Upwind 

0.049 

(0.003) 

0.053 

(0.010) 

0.049 

(0.003) 

for  all  bins  of  the  composite.  As  is  expected  by  parcel 
displacement  theory  for  a  CMASL,  flux  in  all  bins  of  the 
composite  is  down-gradient.  This  property  is  also  found  in 
the  majority  of  updrafts  and  downdrafts  of  Khalsa  and  Greenhut 
(1985) .  Composite  patterns  of  w'r'  resemble,  for  the  most 
part,  those  of  w'w' .  The  largesc  flux  magnitude  is  found  in 
the  middle  center  bin.  Fluxes  along  the  upwind  side  of  the 
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draft  is  slightly  larger  han  those  on  the  downwind  side. 
There  is,  however,  symmetry  along  the  cross-mean-wind  group 
axis. 

The  CD  w'r'  composite  is  presented  in  Table  VI.  As  with 


TABLE  VI 


CD  absolute  humidity  flux  composite  and 
corresponding  standard  deviations. 


w'r' (gs  ^m"^) 
(std.  dev.) 

Left 

Center 

Right 

Downwind 

0.025 

(0.003) 

0.025 

(0.003) 

0.025 

(0.003) 

Middle 

0.023 

(0.004) 

0.034 

0.022 

(0.003) 

Upwind 

0.020 

(0.002) 

0.021 

(0.002) 

0.020 

(0.002) 

the  C'J,  there  is  statistical  confidence  in  the  entire 
composite  and  the  absolute  humidity  flux  is  down-gradient 
throughout  the  composite,  with  a  maximum  at  the  center.  There 
are,  however,  some  differences  from  the  CU  w'r'  composite.  As 
is  expected  from  parcel  displacement  theory  and  the  relative 
magnitudes  of  flux  in  the  w'w'  draft  composites,  a  weaker 
absolute  humidity  flux  composite  is  observed  for  the  CD.  Also, 
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the  downwind  side  of  the  CD  composite  contains  slightly 
stronger  flux  than  that  of  the  upwind  side  but  as  in  the  CU 
composite,  symmetry  exists  along  the  cross-mean-wind  group 
axis. 

3.4  MOMENTUM  FLUX 

There  is  statistical  confidence  for  all  bins  of  the  CU 
planview  along-measn-wind  momentum  flux  composite,  presented  in 
Table  VII.  The  flux  contained  in  all  bins  are  down-gradient, 
given  the  usual  sign  of  the  surface  layer  wind  shear.  This 
same  property  is  seen  in  moist  updrafts  by  Khalsa  and  Greenhut 
(1985).  Wilczak  (1984),  however,  shows  a  region  of  counter¬ 
gradient  momentum  flux  just  inside  the  trailing  edge  of  his 
ensemble  updraft  vertical  cross  section.  The  largest  flux 
magnitude  in  the  CU  composite  is  found  in  the  middle  center 
bin.  Less  significant  asymmetries  exist.  For  example,  w'u' 
is  stronger  along  the  downwind  side  of  the  composite  than 
along  the  upwind  side.  Weakest  flux  is  found  in  the  middle 
left  and  right  bins. 

The  CD  planview  flux  composite  of  w'u'  is  presented  in 
Table  VIII.  Statistical  confidence  exists  in  all  data  with 
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TABLE  VII 

CU  along-mean-wind  momentum  flux  composite  and 
corresponding  standard  deviations. _ 


w'u'  (m^s“^) 

( std .  dev . ) 

Left 

Center 

Right 

Downwind 

-0.133 

(0.008) 

-0.145 

(0.041) 

-0.133 

(0.008) 

Middle 

-0.107 

(0.016) 

-0.219 

-0.099 

(0.018) 

Upwind 

-0.122 

(0.023) 

-0.123 

(0.037) 

-0.122 

(0.023) 

the  exception  of  the  upwind  center  bin,  where  the  standard 
deviation  approaches  the  magnitude  of  the  flux.  As  with  the 
CU,  and  as  seen  in  Khalsa  and  Greenhut  (1985) ,  w'u'  is  down- 
gradient  throughout  the  composite.  As  discussed  below,  this 
finding  differs  from  that  of  Wilczak  (1984)  for  his  ensemble 
downdraft  vertical  cross  section  of  w'u' .  The  region  of 
strongest  flux  within  the  CD  composite  is  found  on  the 
downwind  side,  while  weakest  flux  is  on  the  upwind  side.  The 
largest  magnitude  of  w'u'  is  seen  in  the  downwind  center  bin 
while  the  smallest  is  in  the  upwind  center  bin.  Cross-mean- 
wind  group  axis  symmetry  is  seen. 

As  stated  above,  Wilczak  (1984),  using  finer 
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observational  resolution  than  was  possible  in  this  study, 
shows  a  region  of  counter-gradient  momentum  flux  just  inside 


TABLE  VIII 

CD  along-mean-wind  momentum  flux  composite  and 
_ corresponding  standard  deviations.  _ 


w'u'  (m^s“^) 
(std.  dev.) 

Left 

Center 

Right 

Downwind 

-0.077 

(0.014) 

-0.092 

(0.027) 

-0.077 

(0.014) 

Middle 

-0.065 

(0.015) 

-0.084 

-0.063 

(0.016) 

Upwind 

-0.043 

(0.008) 

-0.028 

(0.022) 

-0.043 

(0.008) 

the  trailing  edge  of  his  ensemble  updrafts  and  downdrafts. 
Given  the  differences  in  horizontal  resolution,  it  is  possible 
that  the  along-mean-wind  asymmetries  in  both  studies  are 


related. 
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4 .  SuBaary 

The  coupling  of  the  atmosphere  to  the  ocean  surface  is 
accomplished  predominately  through  fluxes  induced  by  drafts 
within  the  marine  atmospheric  surface  layer.  Planview 
composite  draft  eddy  correlation  flux  patterns  of  vertical 
velocity,  buoyancy,  absolute  humidity,  and  the  along-mean-wind 
component  of  momentum,  within  the  convective  marine 
atmospheric  surface  layer,  are  derived  in  this  study.  High 
resolution  turbulence  data  are  obtained  from  20  NCAR  Electra 
50  m  mean  sea  level  flight  legs  from  Project  FIRE.  Data 
processing  results  in  a  demeaned,  detrended,  despiked,  high- 
pass  filtered  perturbation  dataset. 

A  conditional  sampling  technique  is  used  to  isolate 
updrafts  and  downdrafts  of  interest  from  the  environment.  The 
technique  uses  horizontal  scale  and  magnitude  thresholds  based 
on  w'  event  criteria.  The  w'  dataset  is  band-pass  filtered  to 
eliminate  the  inertial  subrange  contribution  (0.1  cutoff 
wavelength)  and  the  mesoscale  contribution  (1.5  cutoff 
wavelength)  .  The  minimum  draft  width  is  50  m  while  O.l  m/s  is 
the  w'  magnitude  threshold.  This  technique  results  in  the 
isolation  of  2839  updrafts,  averaging  98  m  in  width,  and  3062 
downdrafts,  averaging  107  ro  in  width. 

After  the  identification  of  drafts  of  interest,  eddy 
correlation  fluxes  are  calculated  for  those  portions  of  the 
processed  dataset  where  drafts  of  interest  are  located.  Mean 
wind  positioning  within  drafts,  draft  symmetry,  and  mean  draft 
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intensity  are  used  to  implement  a  series  of  averaging 
calculations,  grouping  techniques,  and  scaling  operations. 

The  result  of  the  above  data  manipulation  is  a  planview 
composite  updraft  and  downdraft  for  each  flux  discussed.  The 
composites  are  composed  of  flux  data  located  in  horizontal 
three  bin  by  three  bin  arrays.  Each  bin  represents  a  mean- 
wind-relative  sector  of  a  composite  draft:  upwind,  middle,  or 
downwind  by  left,  center,  or  right.  Along-mean-wind  group 
axes  extend  through  the  center  of  composite  drafts  while 
cross-mean-wind  group  axes  extend  through  the  middle  of 
composite  drafts. 

Statistical  confidence  exists  in  the  two  vertical 
velocity  flux  composites.  The  strongest  flux  region  of  the 
updraft  composite  is  located  along  the  along-mean-wind  group 
axis,  with  the  largest  magnitude  flux  found  at  the  center  of 
the  draft.  While  the  largest  magnitude  of  w'w'  within  the 
downdraft  composite  is  also  found  at  the  center  of  the  draft, 
the  overall  composite  is  weaker  and  more  uniform  than  that  of 
the  updraft.  In  the  buoyancy  flux  composites,  only  data  in 
bins  along  the  along-mean-wind  and  cross-mean-wind  group  axes 
are  statistically  significant.  Along  these  axes,  symmetries 
exist  and  the  fluxes  are  down-gradient.  Statistical 
confidence  exists  throughout  both  absolute  humidity  flux 
composites.  As  is  expected  by  parcel  displacement  theory,  all 
fluxes  are  down-gradient.  The  largest  magnitudes  of  flux  are 
found  in  the  middle  center  bins,  and  both  composites  show 
cross-mean-wind  group  axis  symmetry.  The  updraft  and 
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downdraft  along-mean-wind  component  of  momentum  flux 
composites  are  statistically  significant  throughout,  with  the 
exception  of  the  upwind  center  bin  of  the  downdraft  composite. 
All  fluxes  of  both  composites  are  down-gradient.  The  largest 
flux  magnitude  for  the  updraft  is  in  the  middle  center  bin 
while  that  of  the  downdraft  is  located  in  the  downwind  center 
bin.  Along-mean-wind  asymmetry  exists  within  both  composites 
with  downwind  bins  having  greater  magnitudes  of  flux  than 
upwind  bins. 

It  is  hoped  that  this  study  provides  new  useful  insight 
into  the  planview  flux  structure  of  CMASL  updrafts  and 
downdrafts.  The  degree  of  spatial  resolution  in  these 
results,  as  well  as  the  techniques  presented  in  this  study, 
could  be  of  use  in  verification  of  nonlinear  convective 
boundary  layer  models  such  as  that  presented  in  Haack  and 
Shirer  (1992)  .  These  results  are  also  appropriate  for  the 
input  of  forcing  in  time  dependent  two-dimensional  ocean  wave 
models  of  cats  paw-type  features. 

Of  course  much  room  exists  for  improvement  upon  and 
expansion  beyond  this  study.  This  study  should  be  used  as  a 
stepping  stone  for  such  research  initiatives. 
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Abstract 

The  correlation  dimension  d  is  used  commonly  to  quantify  the  chaotic  structure  of 
an  attractor  of  a  smooth  dynamical  system.  The  b^aiidard  algorithm  for  estimating  the  value 
of  d  is  based  on  finding  the  slope  of  the  curve  obtained  by  plotting  In  C(r)  versus  Inr, 
where  C(r)  is  the  correlation  integral  and  r  is  the  distance  between  points  on  the  attractor. 
It  is  argued  here  that  this  algorithm  depends  implicitly  and  sensitively  on  the  assumption 
that  C(r)  is  differentiable,  which  is  an  assumption  apparently  not  applicable  to  even  the 
Lorenz  attractor.  Moreover,  there  is  uncertainty  in  the  value  typically  cited  for  d  because  its 
numerical  approximation  is  often  expressed  as  a  function  of  r  and  because  no  objective 
criterion  is  given  for  locating  the  scaling  region  of  appropriate  values  of  r  to  consider. 
Finally,  a  priori  statistical  error  bounds  on  the  approximate  value  for  d  are  typically  not 
reported. 

In  this  article,  a  rigorous  means  for  obtaining  from  a  single  data  set  an  infinite 
number  of  independent  estimates  for  the  correlation  dimension  d  is  given  that  relies  only  on 
the  three  assumptions  that  C(r)  is  continuous,  that  C(r)  obeys  the  relation  C(r)-r'^  as  r 
approaches  zero,  and  that  a  certain  limiting  moment  exists.  The  most  likely  candidates  for 
d  given  by  this  new  integral  method  are  those  estimates  that  are  relatively  constant  and 
close  to  one  another  in  a  range  of  r.  Moreover,  it  is  demonstrated  that,  when  the  nonlacunar 
relationship  C{r)  =  r‘‘  holds  but  the  empirical  estimate  C,(r)  of  C(r)  is  poor  for  small 
values  of  r,  the  integral  method  converges  more  rapidly  to  the  correct  value  of  d  than  does 
the  standard  slope  method.  Validation  of  the  estimates  for  d  is  obtained  by  making  a  priori 
statistical  error  estimates  that  crudely  bound  the  errors  between  the  approximate  and  actual 
values  of  d.  The  most  likely  value  for  d  is  now  the  one  yielding  the  smallest  error  in  these 
estimates.  Although  these  results  bear  on  the  location  of  the  scaling  region,  the  problem  of 
identifying  that  region  remains  open. 
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A  rigorous  probabilistic  approach  is  used  to  develop  the  expressions  in  this  article. 
For  any  positive  numbers  p  and  p,  the  expected  value  of  {r/pY  over  distances  less  than  or 
equal  to  p  is  denoted  by  E^irj pY  :r<  pj.  The  limiting  moment  M{p)  is  then  defined  to  be 
the  limit  of  this  expected  value  as  p  approaches  zero.  An  infinite  number  of  estimates  for  d 
is  obtained  first  by  varying  the  value  of  p  to  produce  different  values  of  M(p)  and  then  by 
using  the  simple  formula  d  =  pM(p)f{l-M{p)).  A  priori  statistical  estimates  are  given  for 
the  errors  in  the  approximations  of  the  expected  values  E^{r/ pY-f^  p)-  Coarser  estimates 
are  given  for  their  differences  from  M(p),  and  finally  for  the  resulting  approximations  of  d. 
Further  refining  the  identification  of  the  scaling  region  thus  depends  on  improving  the 
estimates  of  the  difference  |  E^(rfpY : r  < p^-M{p)  j . 

In  important  cases,  the  limiting  moment  M{p)  may  not  exist.  But  the  limit  of 
E^{rlpY :r< p^  along  a  suitable  sequence  p  =  p,,p2,p3,...  does  exist,  and  yields  the 
dimension  in  the  same  way.  It  is  a  subject  for  future  research  to  determine  this  sequence. 

Illustrations  using  the  Lorenz  and  Henon  attractors  are  presented.  The  integral 
method  alone  yields  the  following  candidates  for  d:  For  the  Lorenz  attractor,  2.08  or  2.09  is 
indicated;  but  for  the  Henon  attractor,  no  single  value  is  strongly  indicated,  although  the 
values  appear  to  fall  in  the  range  1.21  <</<  1.27.  For  the  Lorenz  attractor,  insufficient 
independent  time  series  points  are  used  to  yield  usable  error  bounds  for  d;  but  for  the  Henon 
attractor,  enough  are  used  to  produce  the  estimate  c?  =  1.25 ±0.15,  which  holds  with  a 
97.5%  probability  and  moment  error  tolerance  of  0.02,  and  the  estimate  «/=  1.25 ±0.1, 
which  holds  with  an  90%  probability  and  moment  error  tolerance  of  0.01.  All  these  values 
fall  within  the  ranges  typically  reported  for  each  attractor. 

A  simple  procedure  for  estimating  the  value  of  d  emerges  from  the  statistical  error 
analysis  and  its  application  to  the  Henon  attractor:  The  optimum  value  for  d  is  the  particular 
candidate  from  the  above  calculation  that  minimizes  the  difference  between  the  bounds  a 
and  A  in  the  inequalities  a  <  C(r)lr‘‘  <  A  that  follow  from  the  strong  version  of  the  relation 
C(r)  ~  r'',  upon  which  the  a  priori  statistical  error  estimates  are  based. 
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Finally,  an  equation  for  the  limit  as  p  approaches  zero  of  E^ir/pY [ln(r/p)]’' :r<p^ 
is  found  that  generalizes  to  the  case  C{r)~r'‘  those  formulas  given  by  Takens  for  the 
nonlacunar  case  C{r)  =  ar‘‘  when  p  =  0  and  y=l  or  y  =  2.  In  addition,  similar  equations 
are  found  for  the  Hentschel  and  Procaccia  generalized  dimensions  with  q  an  integer. 
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1.  Introduction 

Estimates  of  the  correlation  dimension  d  oi  2i  time  series  are  commonly  used  for 
quantifying  the  chaotic  complexity  of  a  physical  system  (e.g.,  Grassberger  and  Procaccia, 
1983a,  b;  Eckmann  and  Ruelle,  1985;  Henderson  and  Wells,  1988;  Baker  and  Gollub,  1990; 
Tsonis,  1992).  However,  owing  to  the  presence  of  both  intrinsic  and  extrinsic  sources  of 
error,  the  merits  of  this  procedure  have  been  debated  considerably  in  the  literature  (e.g., 
Nicolis  and  Nicolis,  1984,  1987;  Grassberger,  1986,  1987;  Smith  et  ai,  1986;  Theiler, 
1986,  1988;  Procaccia,  1988;  Smith,  1988;  Tsonis  and  Eisner,  1988,  1989,  1990;  Nerenberg 
and  Essex,  1990;  Ruelle,  1990;  Essex  and  Nerenberg,  1991). 

Numerous  sources  of  error  inherent  in  the  approximation  of  the  attractor  dimension 
arise  from  such  things  as  time  series  collection  and  preprocessing,  model  reconstruction, 
and  the  specific  dimension  algorithm  used.  Under  time  series  collection,  we  include  the 
frequency,  number,  continuity  and  duration  of  the  measurements  (Nicolis  and  Nicolis, 
1984,  1987;  Grassberger,  1986,  1987;  Fraedrich,  1986;  Smith  et  al,  1986;  Essex  et  ai, 
1987;  Smith,  1988;  Keppene  and  Nicolis,  1989;  Ramsey  and  Yuan,  1989;  Nerenberg  and 
Essex,  1990;  Ruelle,  1990;  Essex,  1991)  and  contamination  by  noise  (Franaszek,  1984; 
Ben-Mizrachi  et  ai,  1984;  Theiler,  1986;  Simm  et  ai,  1987;  Osborne  and  Provenzale, 
1989;  Theiler,  1991);  under  preprocessing,  we  include  low  or  high  pass  filtering  (Theiler, 
1991),  interpolation  and  trend  removal,  and  principal  component  analysis  (Albano  et  ai, 
1988);  under  model  reconstruction,  we  include  the  choice  of  variable  or  variables  (Lorenz, 
1991),  embedding  dimension  (Man6,  1981;  Takens,  1981;  Grassberger  and  Procaccia, 
1983a,  b;  Keppene  and  Nicolis,  1989;  Nerenberg  and  Essex,  1990;  Tsonis  et  ai,  1993),  and 
phase  lag  (Broomhead  and  King,  1986;  Fraser  and  Swinney,  1986;  Albano  et  ai,  1991; 
Thomson  and  Henderson,  1992);  and  under  the  dimension  algorithm  used,  we  include  those 
proposed  by  Kaplan  and  Yorke  (1979),  Russell  et  ai  (1980),  Takens  (1981),  Greenside  et 
ai  (1982),  Grassberger  and  Procaccia  (1983a,  b),  Theiler  (1987),  Henderson  and  Wells 
(1988),  and  Kember  and  Fowler  (1992). 
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In  this  article,  we  focus  on  the  last  general  source,  the  dimension  algorithm  used. 
Thus  it  is  a  principal  objective  of  this  article  to  present,  in  section  2,  new  algorithms  for 
sharpening  the  accuracy  of  the  estimate  for  d.  In  deriving  our  new  algorithms,  we  seek  to 
avoid  those  errors  that  result  from  assuming  that  the  correlation  integral  C(r)  is 
differentiable  when  it  is  not.  For  example,  the  standard  slope  method,  which  uses  the  slopes 
of  the  graph  of  InC(r)  versus  Inr  (e.g.  Grassberger  and  Procaccia,  1983a,  b;  Nicolis  and 
Nicolis,  1984;  Eckmann  and  Ruelle,  1985;  Fraedrich,  1986;  Grassberger,  1986;  Nese,  1987; 
Simm  et  al,  1987;  Henderson  and  Wells,  1988;  Keppene  and  Nicolis,  1989;  Tsonis  and 
Eisner,  1989;  Nerenberg  and  Essex,  1990;  Albano  et  al,  1991;  Lorenz,  1991;  Nerenberg  et 
al,  1991),  rests  implicitly  on  the  hypothesis  that  C(r)  is  differentiable. 

Unfortunately,  careful  numerical  analysis  of  the  Lorenz  attractor  (Lorenz,  1963) 
suggests  that  C(r)  is  not  differentiable.  In  Figure  la,  two  slope  method  estimates  for  d  are 
shown  that  were  each  obtained  using  a  simple  three-point  centered  finite-difference 
approximation.  So  that  problems  introduced  by  employing  model  reconstruction  of  a  single 
series  (e.g.,  Takens,  1981)  were  not  encountered,  both  these  estimates  are  based  on  the  same 
three- variable,  transient-free  time  series  of  20,000  points  separated  by  a  time  interval  of  0.5. 
In  the  first  case,  the  distance  data  were  partitioned  into  500  bins  of  width  0.2,  while  in  the 
second  case,  they  were  partitioned  into  10,000  bins  of  width  0.01.  In  Figure  la,  the 
estimates  for  d  using  data  up  to  bin  distances  r  are  shown;  clearly,  as  the  resolution  is 
increased  by  decreasing  the  bin  width,  the  estimates  for  d  become  noisier,  thereby 
suggesting  that  the  slopes  of  the  correlation  function  are  not  approaching  a  limit;  that  is,  the 
correlation  function  is  not  differentiable.  Similar  results  can  be  inferred  from  those  in 
Cutler  (1991). 


Figure  1 
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Figure  1.  Estimates  of  the  correlation  dimension  d  for  the  Lorenz  attractor.  A  transient- 
free,  three-variable  time  series  containing  20,000  points  is  used  for  the  standard  case 
of  Rayleigh  number  Ra=2S,  Prandtl  number  a=10,  and  domain  shape  parameter 
b  =  Zf3.  In  (a),  the  standard  slope  method  is  compared  when  500  (thick  curve;  bin 

width  0.2)  and  10,000  (thin  curve;  bin  widtii  3.01)  distance  r  bins  are  used;  in  (b)  the 
slope  (thin  curve)  and  integral  (thick  curve;  p  =  l)  methods  are  compared  for  the 
10,000-bin  case;  and  in  (c)  the  integral  method  (p  =  1)  is  compared  using  500  (thick 
curve)  and  10,000  (thin  curve)  bins. 
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If  in  fact  the  function  C(r)  is  not  differentiable,  then  the  standard  slope  procedure 
for  estimating  the  value  of  d  has  no  foundation  and  may  well  break  down.  We  may 
illustrate  such  a  collapse  in  the  slope  procedure  by  considering  the  middle-third  Cantor  set. 
In  this  case,  the  function  C(r)  is  the  well-known  Cantor  function  that  is  almost  everywhere 
differentiable  with  derivative  zero  (Folland,  1984).  For  this  function,  the  slope  procedure, 
especially  if  perfectly  implemented,  would  yield  an  estimate  of  zero  for  the  correlation 
dimension  d  of  the  Cantor  set  rather  than  the  correct  value  of  In2/ln3  (Mandelbrot,  1983; 
Barnsley,  1988;  Cutler,  1991).  The  example  of  the  Cantor  function  presents  deeper 
problems  as  well:  Theiler  (1988)  shows  that  the  Takens  (1985)  formula  fails  to  converge  for 
this  example,  and  below  we  note  that  Theiler's  argument  applies  to  show  that  our  limiting 
moment  M{p)  does  not  exist  for  the  Cantor  function.  We  believe  that  finding  a  rapidly 
converging  algorithm  for  the  Cantor  function  is  an  outstanding  problem  to  solve  in  the 
search  for  improved  dimension  algorithms  (see  also  Ameodo  et  al,  1987). 

In  Figure  lb  we  compare  results  given  by  the  slope  method  and  the  integral  method 
introduced  in  this  article.  Both  estimates  for  d  use  the  same  20,000-point  transient-free  time 
series  from  the  Lorenz  model  that  was  used  to  create  Figure  la,  and  both  use  distance  data 
partitioned  into  10,000  bins.  As  also  illustrated  in  Figure  Ic,  our  method  clearly  does  not 
suffer  from  the  same  noisy  property  as  does  the  slope  method.  Moreover,  as  discussed  in 
section  3,  our  approach  is  not  equivalent  to  merely  averaging  the  slope  estimates  over  some 
interval  of  r,  and  in  fact  gives  a  more  accurate  answer  in  an  idealized  case  when  there  is 
noise  contamination  at  small  values  of  r. 

We  develop  in  section  2  an  integral  method  that  avoids  the  hypothesis  that  C(r)  is 
differentiable  by  integrating  the  product  r'’C(r);  this  method  is  a  generalization  of  a 
dimension  estimator  proposed  by  Theiler  (1988;  1990)  that  in  turn  is  based  on  one  proposed 
by  Takens  (1985).  We  compare  in  section  3  the  accuracy  of  the  integral  method  with  that 
of  the  slope  method  in  a  particular  case  where  C(r)  scales  one  way  {r")  in  an  undersampled 
or  noise-contaminated  region  and  another  way  (r'^)  in  a  well-sampled  region.  For  example. 
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Essex  (1991)  argues  that  undersampling  leads  to  a  value  of  c  that  is  less  than  the  correct 
value  of  the  dimension.  In  contrast,  Ben-Mizrachi  et  al.  (1984)  and  Simm  et  al.  (1987) 
argue  that  the  dimension  c  to  be  expected  from  the  noise-contaminated  region  is  the 
embedding  dimension,  while  the  dimension  from  the  well-sampled  region  is  the  correct  one 
d.  In  our  comparison  in  section  3,  we  find  both  analytically  and  numerically  that  the 
integral  method  converges  much  more  rapidly  to  the  correct  dimension  than  does  the  slope 
method.  In  turn,  Theiler  (1988)  notes  that  the  slope  method  converges  more  rapidly  than 
does  direct  application  of  the  definition  of  the  correlation  dimension. 

Another  advantage  of  the  integral  method  is  that  we  may  interpret  it  in  a  statistical 
way  so  that  we  are  able,  in  section  5,  to  equip  its  dimension  approximants  with  fairly 
natural  a  priori  error  estimates  (c/,  Denker  and  Keller,  1986;  Holzfuss  and  Mayer-Kress, 
1986;  Ramsey  and  Yuan,  1989;  Theiler,  1990;  Judd,  1992;  Kember  and  Fowler,  1992).  We 
apply  in  section  6  the  integral  method  dimension  and  error  estimates  to  the  Henon  attractor 
(Henon,  1976)  and  the  Lorenz  attractor  (Lorenz,  1963). 

Finally,  in  section  4  we  develop  rigorous  generalizations  of  the  Takens  fonriulas  for 
estimating  the  value  of  d  (Takens,  1985)  and  in  section  7  we  derive  a  variant  of  the  integral 
method  for  the  Hentschel-Procaccia  generalized  dimensions  with  q  an  integer  (Hentschel 
and  Procaccia,  1983).  Concluding  remarks  in  section  8  finish  the  article. 


2.  Integral  Method  for  Estimating  d 

We  recall  the  slope  method  for  extracting  the  correlation  dimension  d  from 
knowledge  of  the  correlation  integral  C(r).  First,  we  suppose  that  we  have  the  asymptotic 
relation 

C(r)~r"  (2.1) 

defining  the  dimension  d.  It  follows  immediately  that  the  limit  equation 


'-*0  Inr 


(2.2) 


produces  d.  Unfortunately,  using  (2.2)  is  an  inefficient  way  to  find  d  (Theiler,  1988);  but  in 
the  case  that  C(r)  is  differentiable,  we  may  apply  L’Hopital's  Rule  to  (2.2)  to  obtain  the 
more  efficient  formula  (Theiler,  1986;  Albano  et  al,  1988) 

dlnC(r)  r  dC(x) 

d  =  hm - ^  =  lim - —  (2.3) 

'-*0  dlnr  dr 

Unfortunately  again,  if  C(r)  is  not  known  to  be  differentiable,  then  the  transition  from  (2.2) 
to  (2.3)  is  not  necessarily  valid  and  (2.3)  does  not  make  sense.  Furthermore,  even  if  C(r)  is 
differentiable,  L'H6pital's  Rule  may  not  apply  to  make  the  two  limits  in  (2.3)  equal. 
Accordingly,  we  need  an  alternative  algorithm  in  the  nondifferentiable  case. 

In  this  section  we  develop  such  an  algorithm.  We  begin  by  regarding  the  correlation 
integral  C(r)  as  a  cumulative  probability  distribution  for  the  random  variable  r. 
Consequently,  we  may  define  for  any  positive  number  p  the  pth  normalized  moment  relative 
to  the  condition  r<p  as  the  expected  value  E^{r/  pY  :r  <  pj  of  the  random  variable  {rfpY, 
subject  to  the  condition  r  <  p.  Thus,  we  set 

E{{rlpY  :r<p)  =  —  f  f-1  dCir) 

Cip)jXp) 

Then  we  define  the  pth  limiting  moment  M(p)  of  C(r)  by  setting 

M(p)  =  lim£((r/p)'’ :  r  <  p) 

The  basic  result  of  this  section  is  that,  if  the  limit  in  (2.5)  exists,  then  the  equation 

pM{p) 

~  l-M{p) 

is  valid  for  the  correlation  dimension  d  given  by  (2.2). 


(2.4) 


(2.5) 


(2.6) 
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To  see  that  (2.6)  is  valid,  we  find  it  helpful  to  interpret  the  asymptotic  relation  (2.1) 
in  two  senses,  a  strong  sense  and  a  weak  sense.  We  say  that  (2.1)  holds  in  the  strong  sense 
if  there  exist  positive  numbers  a.  A,  and  Pg  for  which  the  inequalities 

a<C(r)lr‘‘<A  (2.7) 

hold  whenever  the  value  of  r  satisfies 

0<r<Po  (2.8) 

The  numbers  a.  A,  and  Pg  are  parameters  relevant  to  the  statistical  error  estimates  in  section 
5,  where  we  must  therefore  require  the  strong  sense.  We  note  that  the  case  a  =  A  is  the  one 
termed  nonlacunar  by  Theiler  (1988). 

In  contrast,  we  say  that  (2.1)  holds  in  the  weak  sense  if  the  limit  equations 

lim^  =  0  (2.9) 

r-»0  /•'’l 

and 

lim^  =  oo  (2.10) 

r-*0  f‘^2 

are  valid  for  any  two  numbers  and  satisfying 

d^<d<d.  (2.11) 

Of  course,  the  strong  sense  implies  the  weak.  It  is  the  weak  sense  that  suffices  for  the 
derivation  below  of  the  formulas  for  our  new  integral  method. 

We  need  to  make  only  three  assumptions  about  the  function  C(r),  assumptions  we 
believe  to  be  commonly  made  and  reasonable  when  C(r)  is  the  correlation  integral: 
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(i)  We  assume  that  C(r)  is  a  continuous  function  of  r.  That  is,  we  assume 
that  no  particular  value  of  r  is  realized  with  positive  probability. 

(ii)  We  assume  that  the  asymptotic  equation  C(r)~r‘'  holds  in  the  weak 
sense. 

(iii)  For  some  value  of  p*0  such  that  d+p>\,  we  assume  that  the  limiting 
moment  M{p)  of  (2.5)  exists. 


The  argument  that  these  three  assumptions  yield  (2.6)  is  embodied  in  an  easy 
calculation.  First  we  integrate  by  parts  to  obtain 


j%'’dC{r)  =  p^C{p)-pj^^r'’  ^C(r)dr 


(2.12) 


We  combine  (2.12)  with  the  definition  (2.4)  to  find 

,  .  {'r^-'C{r)dr 

eUvIpY  ;  r  <  p)  =  1  -  p— — — - 

y  p^C(p) 

Therefore,  because  the  limit  M{p)  exists  by  assumption  (iii),  we  obtain 

rr'’-‘C(r)t/r 

A/(p)  =  1  -  plim  - 

pPCip) 


(2.13) 


(2.14) 


in  which  it  follows  that  the  limit  on  the  right  side  exists.  Next  we  verify  that  we  have 

j%^’C(r)ifr (2.15) 


in  the  weak  sense,  from  which  it  follows  that 


In  r 

p+d  =  \m—^ - 

Inp 


(2.16) 


Because  we  have  assumed  that  d+p>l,  we  see  that  the  function  r'’”'C(r)  of  r  is 
continuous,  and  then  that  Yr'^^Cir}dr  is  a  differentiable  function  of  p.  Consequently,  we 
may  apply  L'Hdpital's  Rule  to  obtain 
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=  (2.17) 

P-o  f%^ic{r)dr 

Jo 

in  which  the  limit  exists  by  (2.14)  and  assumption  (iii).  Finally,  we  obtain  from  (2.14)  and 
(2.17) 


M{p)  =  l-p 


1 

p+d 


d 

p+d 


(2.18) 


from  which  (2.6)  follows. 

We  note  that  in  the  case  p  =  0,  the  same  Stieltjes  integration  by  parts  argument  leads 
to  the  Takens  estimator  (Takens,  1985) 


f’-dC(r)  =  -E(lHr/p)]-.rSp)  =  ^  (2,19) 

C(p)  J  0  r  d 

which  is  essentially  the  content  of  (2.15)  and  (2.16)  in  Theiler  (1988)  or  of  (22)  and  (23)  in 
Theiler  (1990). 


Remark  1.  There  is  an  important  subtlety  at  this  point:  L’Hopital’s  Rule  does  not 
imply  that  if  the  limit  in  (2. 16)  exists,  then  so  does  M(p).  In  fact,  the  limit  in  (2.16) 
exists  whenever  C(r)  ~ /  holds,  even  in  the  case  p  =  0,  but,  as  Theiler  (1988)  notes, 
the  corresponding  moment  Af(0, 1)  defined  in  (4.19)  does  not  exist  for  C(r)  the 
Cantor  function;  accordingly,  the  Takens  estimator  (2.19)  fails  in  that  case.  In  fact, 
Theiler's  argument  applies  to  show  that  also  the  limit  M(p)  fails  to  exist  for  C(r) 
the  Cantor  function. 


Remark  2.  Even  if  the  limit  in  (2.17)  fails  to  exist,  the  Cauchy  Mean  Value 
Theorem  (used  to  prove  L'HSpital's  Rule;  James,  1967)  implies  that  there  exists  a 
sequence  p, ,  Pj ,  p,, . . .  of  values  of  p  converging  to  zero,  so  that  the  limit 
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M,  ip)  =  lira  £((r/p*  Y:r<p,) 

k-*oa  '  ' 

exists.  Then  the  correlation  dimension  is  given  by 


pMAp) 

(l-MsiP)) 


(2.20) 


(2.21) 


3.  Comparison  of  Slope  and  Integral  Methods  for  Estimating  d 

To  form  some  idea  of  the  relative  accuracies  of  the  slope  and  integral  methods,  we 
consider  a  case  in  which  the  nonlacunar  correlation  integral  C(r)  is  approximated  by  the 
empirical  correlation  function 


for  0<r<5 
jr''  for  8<r<p 


(3.1) 


where  5  is  a  positive  number;  this  function  is  given  by  the  solid  curve  in  Figure  2a. 

We  take  (3.1)  as  a  paradigm  for  the  common  situation  in  which  small  amplitude 
noise  distorts  the  empincal  determination  of  the  correlation  integral  C(r)  for  small  values  of 
r.  In  the  idealized  situation  represented  by  (3.1),  we  postulate  that  we  recover  C(r) 
perfectly  for  values  of  r  in  the  range  5<r<p,  while  small  amplitude  noise  leads  to  a 
growth  exponent  c  different  from  the  actual  value  d  in  the  range  0<r^5.  Ben-Mizrachi  et 
al.  (1984),  Simm  et  al.  (1987)  and  Theiler  (1986,  1991)  have  argued  that  an  empirically 
determined  correlation  function  may  be  expected  to  scale  as  r\  with  c  the  embedding 
dimension,  for  small  values  of  r  where  noise  contamination  dominates,  and  as  r'',  with  d  the 
correlation  dimension,  for  larger  values  of  r.  Thus,  following  these  studies,  we  select  c>d 
as  the  most  realistic  case  for  our  example.  The  choice  c  =  2d+l  (Mand,  1981;  Takens, 
1981)  guides  our  illustration  of  C,(r)  in  Figure  2a  with  d  =  2  and  c  =  5. 


Figure  2 


power 
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Figure  2.  Comparison  of  slope  and  integral  methods  for  the  nonlacunar  case  C(r)  =  /  that 
is  approximated  by  the  empirical  correlation  function  C^(r)  in  (3.1),  with  d  =  2  and 
c  =  5.  Data  are  well  sampled  in  the  working  interval  [Spp],  but  are  noise- 
contaminated  in  the  interval  [5i,5].  In  (a),  C,(r)  (solid  curve)  and  the  curve  for  the 
best-f/  rlope  m(5,,p),  (3.4),  (dashed  curve)  are  shown  for  the  case  5]  =0.6,  5  =  0.7, 
and  p  =  0.9.  In  (b)  and  (c),  the  approximate  dimensions  given  by  the  slope  m(5i,p), 
(3.4),  (dashed  curve)  and  integral  d{p,5^,p),  (3.9),  (solid  curve)  methods  are 
compared.  In  (b),  the  case  5,  =0.04,  5  =  0.2,  and  p  =  ^  is  illustrated  as  a  function  of 
p.  In  (c),  the  case  5,  =5^  p  =  j,  and  p  =  0.9  is  shown  as  a  function  of  5.  In  both 
these  cases,  the  integral  method  converges  faster  and  is  less  sensitive  to  noise 
contamination  than  is  the  slope  method. 
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Of  course,  in  actual  cases,  there  will  be  an  intermediate  region  in  which  the 
empirical  correlation  function  varies  from  being  completely  inaccurate  for  small  values  of  r, 
to  being  essentially  accurate  for  larger  values  of  r.  In  our  example,  we  eliminate  this 
intermediate  region  for  simplicity  of  calculation:  The  noise-dominated  region  gives  way 
abruptly  to  the  well-sampled  region  at  r  =  5.  Again,  in  actual  cases,  there  seems  to  be  no 
algorithm— or  even  advice— for  deciding  where  the  noise-dominated  region  ends  and  the 
well-sampled  region  begins.  We  incorporate  this  fact  into  our  example  below  by  not 
allowing  knowledge  of  the  value  of  5  to  enter  into  our  line  fitting  or  integrating  procedures; 
we  use  the  value  of  5  only  to  evaluate  the  result  of  these  procedures.  However,  we  consider 
that  the  working  interval  [S^p]  of  r  used  to  approximate  C(r)  includes  both  well-sampled 
and  noise-dominated  regions.  Thus  we  hypothesize  that  the  inequalities 

0<5,<5<p<l  (3.2) 

are  valid;  the  left-most  inequality  represents  the  fact  that  in  any  sampling  of  a  time  series, 
there  is  a  smallest  positive  distance  5,.  We  expect,  of  course,  that  the  accuracy  of  each 
calculation  will  improve  as  the  length  5  of  the  error-contaminated  interval  [0, 5]  decreases 
in  comparison  with  the  length  of  the  working  interval  [5,,p].  Slightly  more  roughly,  we 
expect  the  accuracy  to  improve  as  5/p  decreases  in  magnitude.  Thus,  we  compare  the  slope 
and  integral  methods  by  discovering  the  rate  at  which  this  accuracy  increases  as  the 
magnitude  of  5/p  decreases. 

To  implement  the  «lope  method  in  the  standard  way,  we  calculate  the  slope  m  of  the 
line  y  =  m\nr  +  b  that  best  fits  the  curve  y  =  lnC,{r)  between  r  =  5,  and  r  =  p  in  y-lnr 
space  (Holzfuss  and  Mayer- Kress,  1986).  If  we  measure  the  closeness  of  fit  by  means  of 
the  mean-square  integral 


l{m,b)=  r  (m]nr  +  b-lnC^(r)f  dr 

JSt 


(3.3) 
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then  we  find  that  the  value  m  =  m{5^,p),  belonging  to  the  minimum  of  is  given, 

with  the  help  of  the  symbolic  manipulator  Derive,  by 


{ S  5^ 

m(5,,p)-</  =  O  — In—  (3.6) 

VP  P) 

Of  course,  when  5,  =5  we  obtain  from  (3.4)  the  expected  equality  m{5^,p)-d  =  0. 
However,  so  long  as  our  fitting  interval  [5, ,  p]  is  contaminated,  even  slightly,  by  the  noise- 
dominated  region,  the  difference  between  the  approximate  dimension  m(5,,p),  obtained 
from  the  slope  method,  and  the  actual  dimension  d  approaches  zero  only  at  the  rate  of 
{5/p)\n{5lp)  as  the  ratio  5/p  approaches  zero. 

The  rate  of  convergence  as  either  p  or  5  varies  is  illustrated  in  Figures  2b,  c,  in 
which  the  approximate  dimension  m(5,,p)  is  denoted  by  the  dashed  curves.  In  Figure  2b,  5 
and  5,  are  held  fixed  and  the  approach  of  m{5^,p)  to  d  shown  as  the  value  of  p  increases. 
In  Figure  2c,  p  is  held  fixed  with  5^  =  5^,  and  the  approach  of  w(5,,p)  to  d  illustrated  as 
5  — >  0 .  In  both  cases,  the  convergence  is  seen  to  be  quite  slow,  meaning  that  even  a  small 
interval  of  r  in  which  there  is  noise  contamination  introduces  a  considerable  error  in  the 
slope  dimension  estimate. 
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In  contrast,  the  difference  between  the  approximate  dimension  obtained  from  the 
integral  method  and  the  actual  dimension  d  approaches  zero  much  faster,  at  least  when  the 
sum  p+d  is  greater  than  1,  as  required  by  assumption  (iii)  in  section  2.  More  specifically, 
we  use  as  the  interval  of  integration  the  same  working  interval  as  above,  [5,,p],  and  we  use 
the  formula 


1 

M(p,S„p)  =  ~ -  -  dC,{r)  (3.7) 

to  define  the  corresponding  approximate  pth  moment  M(p,  5,,p).  Then  we  use  (2.6)  to 
define  the  corresponding  approximate  dimension  d{p,8^,p)  &s 


d{p,5^,p)-d 


(3.10) 


Thus,  as  the  ratio  5/p  approaches  zero,  we  find  that  the  integral  approximate  dimension 
d{p,5^,p)  approaches  the  actual  dimension  d  much  faster  than  the  slope  approximate 
dimension  m{5^,p)  (3.4)  does  when  p+d>l  holds.  This  rapid  approach  is  similar  in  spirit 
to  that  of  Theiler  (1988)  who  noted  that  his  integral  method,  which  is  similar  to  ours  for 
p  =  0 ,  converges  more  rapidly  than  the  slope  method  to  the  actual  dimension  as  p  ->  0 .  For 
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p^O,  the  rapid  approach  is  illustrated  by  the  solid  curves  in  Figures  2b,  c  for  the 
conservative  case  p  =  -j.  Finally,  we  note  from  Figure  2b  that  if  we  interpret  the  interval 
[5i,  5]  as  an  undersampled  region,  then  we  obtain  an  underestimate  for  d  in  that  region,  in 
agreement  with  Essex  (1991). 

In  contrast  to  the  slope  method,  we  note  that  the  integral  method  yields  an  error  in 
the  approximate  dimension  d{p,8^,p)  even  when  there  is  perfect  sampling  in  the  range 
[5i,p]  because  the  integral  (3.7)  upon  which  the  dimension  estimate  is  based  really  should 
be  calculated  over  the  entire  range  [0,  p].  Thus  the  limiting  value  of  d{p,  5, ,  p)  when  5=5, 
and  c  =  d  is  not  d  but 


d(p,5i,p)  =  d  +  {p+d)- 


VP. 


p+d 


Ip. 


(3.11) 


However,  this  result  does  not  affect  the  conclusion  given  by  (3.10)  concerning  the 
convergence  rate  as  the  noise-dominated  region  shrinks  in  size. 


4.  Extension  of  Moment  Formulas  to  Those  of  Takens 

In  this  section  we  derive  the  family  of  formulas  given  by 

lim  £((r/p)"  [ln(r/p)]’' :  r  <  p)  =  (- 1)’'  (4. 1 ) 

(p+d) 

Takens  (1985)  obtained  the  two  estimators 

Urn  £([ln(r/p)];  r  <  p)  =  -3  (4.2) 

P-.0  (J 


and 
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Um£([ln(r/p)f  :r<p)  =  -^  (4.3) 

for  the  special  case  of  a  nonlacunar  correlation  function  C(r)  that  scales  perfectly  with  r; 
that  is,  C{r)  =  ar^ .  Here  we  extend  (4.2)  and  (4.3)  to  all  of  (4.1)  by  deriving  (4.1)  for  the 
general  case  C(r)  -  r**,  where  we  interpret  this  asymptotic  relation  in  the  weak  sense  defined 
in  section  2. 

We  begin  our  derivation  by  noting  that  we  may  interpret  some  of  the  calculations  of 
section  2  as  describing  the  effect  of  rescaling  on  suitable  measures  on  [0, 1].  We  suppose 
that  p  is  a  nonatomic  probability  measure  on  [0,1]  with  the  property  that  every 
neighborhood  of  0  in  [0, 1]  has  positive  p-measure.  A  particular  example  is  the  measure  Vj 
defined  by  setting 

(4.4) 

where  X  is  a  Borel  subset  of  [0, 1].  Given  such  a  measure  p,  for  any  value  of  p  satisfying 

0<p<l  (4.5) 

we  may  define  another  such  measure  Pp  by  setting 


Hp(X) 


p[0,p] 


Using  this  notation,  we  may  formulate  an  easy  theorem: 


(4.6) 


Theorem  If  p  is  a  nonatomic  probability  measure  satisfying  the  condition 

C(p)  =  p[0,p]-p"  (4.7) 

for  small  values  of  p  and  possessing  limiting  moment  M(p)  for  p  =  0, 1, 2, . . .,  then 
the  limit  equation 

'“/p  =  (4.8) 

is  valid  in  the  weak  topology  for  probability  measures. 
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PROOF  We  calculate  the  moment  Is"  dn  As)  of  the  measure  /r  .  We  have 


in  which  s-rjp  and  dp^is)  =  dCir)ICip).  Consequently,  we  may  use  (2.4),  (2.5) 


and  (2.18)  to  obtain 


limfs'’4i  (5)  = 

o-*oJ  ^ 


d+  p 


However,  for  our  standard  measure  Vj,  we  see  easily  that 


js^dvAs)  = 


d+  p 


for  any  value  of  p  >  0.  In  particular,  we  have 


(4.10) 


(4.11) 


limjs'’i/p^(s)  =  js'’tfi>^(5) 


(4.12) 


Consequently  we  find,  for  any  polynomial  P{s),  that 


lim  J/»(j)rfpp(j)  =  I  P{s)dvAs) 


(4.13) 


and  therefore  for  any  continuous  function  f{s),  that 


(4.14) 


Then,  by  definition  of  the  topology  on  the  space  of  Borel  measures  on  [0, 1],  we 


limpp  =  Vj 


(4.15) 
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As  a  corollary  of  this  theorem  (which  actually  is  the  next-to-last  step  in  its  prooO. 
we  obtain  the  limit  equation 

nmJ/(s)<i^p(j)  =  J/(5)rft>^(i)  =  (4.16) 

that  is  valid  for  any  continuous  function  f(s).  Noting  that 

£(/(r/p):  r  Sp)  =  J  f{s)dlip{s)  (4. 17) 

holds,  we  obtain 

lim£(/(r/p):r<p)  =  if|^  f{s)s‘^~^ds  (4.18) 

Now,  with  /(s)  =  ^'’(ln5)’'  and  p+d>l,  we  finally  obtain  the  extension  of  the  Takens 
(1985)  formula  for  a  generalized  limiting  moment: 

A/(p,  y)  =  lm£((r/p)'’[ln(r/p)]’':r  <p)  =  </£5'’[ln5]’' (-1)’'—^^  (4.19) 

This  equation  yields  our  formula  (2.18)  in  the  case  y  =  0  and  includes  the  Takens  (1985) 
formulas  as  the  cases  p  =  0,  y  =  1  and  p  =  0,  y  =  2,  which  are  given  by  (4.2)  and  (4.3). 

In  practice,  as  outlined  in  Takens  (1985),  (4.2)  gives  us  a  way  to  estimate  the  value 
of  J  as  a  mean  of  a  random  sample,  while  (4.3)  gives  us  a  way  to  estimate  the  associated 
variance  by  means  of  the  Kolmogorov  Inequality.  In  the  following  section  we  make  use  of 
the  parallel  observation  for  E^j^rj pY :  r  <  p j  to  arrive  at  statistical  error  estimates. 
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5.  A  Statistical  Relative  Error  Bound  for  d 

Takens  (1985)  notes  that  the  special  case  p  =  Q,  y=l  of  (4.19)  (essentially) 
expresses  the  mean  of  a  random  variable  with  variance  (essentially)  expressed  by  the  special 
case  p  =  0,  7=2.  He  also  observes  that  it  is  therefore  possible  to  incorporate  statistical 
methods  into  an  algorithm  for  estimating  the  value  of  the  correlation  dimension  d.  More 
recently,  Theiler  (1990)  and  Kember  and  Fowler  (1992)  use  such  methods  to  determine 
statistical  error  bounds  on  d. 

Here  we  seek  a  priori  statistical  error  bounds  for  estimates  of  the  correlation 
dimension  d  given  by  the  integral  method.  We  first  apply  the  conventional  use  of  mean  and 
variance  to  obtain  quick  but  probably  unreliable  error  bounds.  Then  we  derive  more 
reliable  error  bounds  with  the  objective  of  identifying  explicitly— and  thereby  making 
computationally  accessible— the  underlying  assumptions  necessary  to  derive  such  bounds. 

To  apply  our  knowledge  of  mean  and  variance  to  find  conventional  error  bounds  for 
an  estimate  of  d,  we  envision  a  situation  in  which  we  may  make  a  large  number  n  of 
observations  of  the  random  variable  (rfpf  for  a  fixed  value  of  p.  We  assume  that  the  value 
of  p  is  so  small  that  the  value  of  E{{rfpY:r<p^  is  essentially  indistinguishable  from  that 
of  M{p)  so  that  we  may  interchange  M(p)  and  E^irf pY :  r  <  p^  in  our  calculations.  Then 
we  may  regard  M{p)  as  the  expected  value  of  our  random  variable  {rjpY ,  and  so  regard 
V{p)  =  M(2p)- M{pf  as  its  variance.  We  next  assume  that  our  n  observations  are 
independent  so  that  the  mean  of  these  observations  has  a  variance  of  V{p)/n.  Using  the 
Central  Limit  Theorem  (Feller,  1966),  we  know  that  when  the  value  of  n  is  large  enough, 
the  mean  of  n  independent  observations  is  essentially  normally  distributed;  we  therefore 
assume  that  in  our  case,  the  value  of  n  is  this  large. 

Accordingly,  the  95%  confidence  interval  J  for  M{p)  is  given  by 

J  =  ^Mip)-2^{p)/n,M{p)  +  2ylV{p)/n^  (5.1) 

Using  (2.18)  for  M{p),  we  may  rewrite  (5.1)  as 
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J=  -^fl-2  I  ,  +  2  I  ----  -  (5.2) 

p+d[  }ln{2p+d)  j  p+d[  \n{2p+d)j 

Thus  we  conclude  that,  with  p  sufficiently  small  and  n  sufficiently  large,  the  probability  is 
0.95  that  the  mean  of  n  independent  observations  of  (r/pY  is  within  a  relative  error  of 
±2^/[nd{2p+d)]  of  the  actual  value  of  M{p).  Then,  upon  applying  (2.6)  to  our 
estimate  (5.2)  of  M(p),  we  find,  with  probability  0.95,  the  value  of  the  correlation 
dimension  d  within  a  relative  error  of  ±2^j{p+df  l\nd{2p+d)^ .  For  example,  if  we  have 
d  =  2,  p  =  3  and  /i  =  10*,  then  we  should  recover  d  within  a  relative  error  of  ±0.0025,  so 
that  if  =  2± 0.005,  with  probability  approximately  0.95.  Of  course,  such  high  precision  is 
not  attainable  in  practice  (e.g..  Fig.  Ic;  Grassberger  and  Procaccia,  1983a,  b;  Abraham  et 
ai,  1986;  Nese  etai,  1987;  and  Simm  etai,  1987). 

The  failure  of  conventional  means  to  give  realistic  error  bounds  for  an  estimate  of 
the  correlation  dimension  suggests  strongly  that  the  values  of  p  normally  used  are  not 
sufficiently  small  to  allow  us  to  interchange  such  quantities  as  M{p)  and  E{{rlpY  :r<pY 
that  the  observations  of  quantities  such  as  {r/pY  are  not  sufficiently  independent  to  allow  us 
to  simply  divide  the  variance  of  the  mean  by  the  number  n  of  observations,  or  that  the 
number  of  independent  observations  is  not  sufficiently  large  to  justify  our  appeal  to  the 
Central  Limit  Theorem.  Of  these  three  difficulties,  perhaps  the  most  serious  is  the  first,  that 
p  is  not  small  enough,  as  emphasized  by  Caswell  and  Yorke  (1986). 

In  view  of  these  difficulties,  we  use  much  more  elementary  and  direct,  but  laborious, 
methods  in  the  remainder  of  this  section  to  obtain  realistic  a  priori  statistical  relative  error 
bounds  for  the  estimates  of  d  given  by  the  integral  method.  That  is,  for  d^  an 
approximation  obtained  in  a  canonical  way  from  a  sample  of  the  attractor,  we  find  a 
formula  bounding  below  the  probability  that  the  inequality 


(5.3) 
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holds.  In  our  derivation  of  that  probability  formula,  we  find  it  convenient  to  alter  the  usual 
definition  of  the  function  C(r)  to  an  essentially  equivalent  definition.  In  the  usual 
definition,  we  consider  a  generic  trajectory  of  points  on  the  attractor.  For  k 

any  positive  integer,  we  calculate  the  distances  for  /<  j<k,  and  we  let  these 

distances  be  r, ,  Tj,  . . .,  r,,  in  order  of  calculation,  with 


n  =  — - - 

2 


(5.4) 


We  set 


=  j^card|(/,7)||jr. -x^|<r  and  i<j<k^  (5.5) 

where  card(A')  denotes  the  number  of  points  in  set  X.  Because  the  sequence  jr, ,  x,,  j:,,  . . .  is 
generic,  we  may  assume  that  the  limit 

C(r)=  limCj(r)  (5.6) 

k-^oo 

exists  and  so  defines  the  correlation  function  C(r).  In  our  altered  definition,  we  consider  a 
generic  pair  of  trajectories  x,,  Jtj,...  and  y,,y2,...  and  we  define  the  time  series  r,,r,,...  of 
distances  by  setting  /;  =  |jr_  -y,  |.  Then  we  redefine 

C,(r)  =  — card{/|/;  <r  and  /<  j<n}  (5.7) 

n 

and 

C(r)=  limC,(r)  (5.8) 

l•-♦0O 

This  view  is  more  convenient  because  we  need  not  renumber  the  series  r,,r2,...,r,  for 
n  =  k{k-\)l2  with  each  successive  choice  of  k.  More  important  is  the  fact  that  thi.s 
interpretation  avoids  the  implicit  dependence  via  triangle  inequalities  of  the  distances 
I  jr,  -xj  (Theiler,  1990);  this  interpretation  also  makes  it  possible  below  to  measure  the  near 
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independence  of  the  distances  r,  by  means  of  autocorrelation.  However,  the  two  views  are 
equivalent  and  our  results  below  using  the  altered  definition  hold  also  for  the  original  one. 
For  each  value  of  p,  we  define  the  number  N{p,n)  of  distances  less  than  p  by 

setting 


N{p,n)  =  '^XirJp) 
1=1 

where  is  the  characteristic  function  defined  by 


fl  if  141  <1 
1 0  otherwise 


Then  we  define  the  approximate  pth  moment  M{p,p,n)  as 


^^(AP.»)  =  -7:  * Xin/P) 
N(p,n)ti\p) 


The  approximate  dimension  is  then  given  by  definition  as  (cf.  (2.6)) 


(5.9) 


(5.10) 


(5.11) 


=  dip,p,n)  = 


pM{p,p,n) 

\-M(p,p,n) 


(5.12) 


With  these  definitions,  we  derive  our  error  estimate  in  three  major  steps.  We  need 
to  assume  that  C{r)~r‘‘  holds  in  the  strong  sense  with  parameters  p^,,  a  and  A  (see  section 
2).  In  our  first  step,  we  obtain  the  coarse  absolute  bound 


E({rfpy:r<p)-M{p) 


<  p 

A  a 

p  +  d 

a  A 

(5.13) 


As  we  see  below,  this  is  not  a  difficult  bound  to  obtain.  In  our  second  step,  we  obtain  a 
simple  Kolmogorov  error  bound 

Prob(|A/(p,p,n)-£((r/p)'’:r  <p)|  ^  t)S:  1-9  (5.14) 
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where  <p  depends  on  n,  r,  a  and  A.  In  (5.14),  the  symbol  ^  means  "approximately  greater 
than"  in  the  sense  explained  below.  Finally,  we  combine  (5.13)  and  (5.14)  with  (5.12)  and 
(2.6)  to  obtain  the  following  theorem: 


StatisticaJ  Relative  Error  Bound  Suppose  that  we  have  positive  numbers  A,  a  and 
P(,  so  that  the  inequality 

a/ ^C{r)<Ar‘‘  (5.15) 

holds  for  0  <  r  <  Po .  In  addition,  suppose  that  0<T<^,0<<p,  0<p<Po  and 


n  >  nip,  T,  p)  = 


40  {2A^+aA) 
9  a^pWp 


hold.  Then  the  inequality 


^e(Ai,x,p,dJ  ^\-p 


holds,  where  we  set 


(5.16) 


(5.17) 


and 


e(Ai,  x,p,dj  = 


A^p+xjp+dJ  p+d^ 
d„-ip+dJ{A^  +  x)  p 


A 

a 


a 

~A 


(5.18) 


(5.19) 


More  precisely,  (5.17)  should  read 


Prob 


d-d„ 


<e(A,,T.p,</J  >l-<p- 


80  Tj 

9  p^^x^a^ 


(5.20) 
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where  r]>0  is  a  measure  of  the  interdependence  of  the  random  variables 
discussed  below;  the  value  of  t]  may  be  made  as  small  as  desired  by  appropriately  lagging 
the  two  time  series  and  yi.y^’  •  •.v«- 

As  outlined  above,  our  first  step  in  deriving  (5.17)  is  to  verify  (5.13).  We  may  use 
the  following  simple  argument  to  estimate  the  limiting  moment  M(p)  with  some  degree  of 
accuracy.  We  recall  that  the  moment  E^{rlpY :r<pj  is  given  by  the  Stieltjes  integral  (2.4); 

£((r/p)':rip)  =  ^j%'rfC(r)  (5.21) 

Then  integration  by  parts  leads  us  to 


£({r/p)' :  r  <  p)  =  I  -  ^ f  r-  C(r)Jr  (5.22) 

Finally,  for  0  <  p  <  Pg ,  we  may  use  (2.7)  to  conclude  that  we  have 


a  p+d 


^E({r/pY:r<p)<\-^ 


A  p  +  d 


(5.23) 


Thus  we  see  that  for  0<p<p„,  both  the  moment  E[i^r  j  pY r  <  p^  and  its  limit 
A/(p)  as  p— >0  must  lie  in  the  interval  [l -(y4/o)(p/(p  +  t/)),  1 -{a/.4)(p/{p  +  c/))j. 
Consequently,  we  have  the  bound  (5.13)  for  the  error  j£((r/p)'’:r  <  pj-M(p)  .  In  general, 
this  is  a  very  crude  bound;  however,  this  bound  is  not  crude  when  the  values  of  A  and  a 


nearly  agree,  as  in  the  case  of  the  Henon  map  discussed  in  the  next  section. 

Our  first  step  completed,  we  turn  to  the  second  one  of  deriving  (5.14).  As  noted  at 
the  beginning  of  this  section,  we  study  the  distances  between  corresponding  points  on  a  pair 
of  time  series  in  the  attractor.  More  specifically,  we  consider  a  dynamical  sy.stem  given  by 
an  invertible  C‘  diffeomorphism  F:  -+  with  a  chaotic  compact  attractor  that  is  the 
closed  support  of  a  probability  measure  p,  invariant  and  ergodic  with  respect  to  F.  In 
addition,  we  suppo.se  that  the  measure  pxp  on  R^xR'^  is  strongly  mixing  and  therefore 
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ergodic  with  respect  to  the  map  FxF.  We  let  F'  =FoF‘  '  and  F^  =F.  Then  by 
definition,  for  x/i -almost  all  initial  point  pairs  (jfo.yo)*  series  of  distances 

;;=|f;Co-Fyo|  (5.24) 

may  be  used  to  compute  the  normalized  moment  E^{r/ p)'’ :  r  <  p j  as  the  limiting  average 


£((r/p)':r<p)=l|m 


1 


N{p,n)‘ 


.p) 


X{rjp) 


(5.25) 


where  x(0  is  the  characteristic  function  defined  by  (5.10)  and  N{p,n)  is  the  number  of 
distances  r.  such  that  r,  <  p  and  i  <  n,  defined  by  (5.9). 

For  simplicity,  we  assume  that  each  new  value  of  {rJpY  xirjp)  is  an  independent 
estimate  of  the  mean  in 


Mo(p,p,n)  =  -^-Xf^'j  Xin/P) 

Using  our  ergodic  hypothesis,  we  note  that  we  may  write 


(5.26) 


E({f/PY  Xif/P))  =  ’  P-  (5.27) 

As  we  show  next,  this  expectation  is  closely  related  to  E^{r/ pY :  r  <  p^,  the  one  we  seek. 
Using  our  ergodic  hypothesis  again,  we  see  that  we  have 


C(p)  =  lim 


N{p,n) 


(5.28) 


thereby  allowing  us  to  take  the  limit  in  (5.25)  with  the  result 


E({r/pY  :r<p)  =  —  E((r/pY  x('-/p)) 


(5.29) 
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To  see  that  in  a  fairly  general  sense  we  may  assume  that  each  new  value  of 
i^ilpY  Xi’^Jp)  be  regarded  as  an  independent  estimate,  we  note  that  the  mixing 
hypothesis  implies  that  the  formula 

Jim  J  J  fix,  y)f{F'”x,  F”y)dn{x)diiiy)  =  (J  j  fix,  y)dnix)dmy)f  (5.30) 

holds  for  any  (integrable  and  square  integrable)  function  fix,y).  We  apply  (5.30)  to  the 
two  choices  fix,y)  =  i\x-y\lpyx^x-y\lp)  and  /(x,y)  =  ;f(|jc-y|/p)  to  find  a  value  of  m 
making  the  difference  of  the  resulting  integrals  smaller  than  a  preset  positive  number  t]  for 
all  values  of  p  smaller  than  son.e  p^.  Then  we  replace  F"”  with  F  in  the  result  so  that  we 
have 

(J  \fix,y)dpix)dpiy))  -  TJ  ^  j  \fix,y)fiFx,  Fy)dpix)dpiy)  <  (j  j  fix,  y)dpix)dpi\ ))“  +  r] 

(5.31) 

for  our  two  choices  of  fix,  y).  Then  (5.31 )  expresses  the  near  independence  we  seek. 

The  conditional  variance  Var(p,  p)  of  the  random  variable  {rjpY  x{^lp)>  subject  to 
the  constraint  that  r  <  p,  is  given  by 

Var(p,  p)  =  E[irlpY'’xirlp))  -  [E(irlpY  x{r Ip))]  (5.32) 

Then  we  may  use  the  Kolmogorov  Inequality  (Feller,  1957)  to  bound  the  probability  that 
the  approximate  mean  Mo(p,p,n),  given  by  (5.26),  differs  from  the  actual  value 
E[{rlpY X{r Ip))  by  a  magnitude  greater  than  p‘‘*^,  wht  3  5>0  is  to  be  determined.  That 
is,  we  have 

Prob(|M,  (p,  p.  n)  -  E((rlpY  xirlp))]  >  p"*)  £  ^  (5.33) 

where  rj  is  the  deviation  from  independence  associated  with  (5.31)  for  our  two  choices 
fix,y)  =  i\x-y\lpY X^x-y\lp)  and  fix,y)  =  x{\x-y\lp)  ■  We  may  rewrite  (5.33)  in  the 
approximate  form 
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Prob(|  W.  (p,  p.  n)  -  £((r/p)'  x(rlp))\  2  <5-34) 

In  the  same  way,  we  note  that  we  have,  by  definition, 

C(p)  =  E{xir/p))  (5.35) 

Then  we  see  that,  because 

X{rl  pf  =  X{rl  p)  (5.36) 

holds,  we  must  have 

C{p)-C{pf  =  Var(x{r/p))  (5.37) 

Then  we  may  use  the  Kolmogorov  Inequality  again  to  obtain,  approximately 

Prob(l Nip, n)/n - C{p)|  > p^*^ ) <  (5.38) 

P  n 

We  may  combine  (5.33)  and  (5.38)  to  obtain  the  estimate 

Prob(|  Mo  ip,  p,  n)  -  E{ir/pY xir/p)]\ <  p'*^  and  |  Nip,  n)/n  -  C(p)|  <  p‘'*^) 

^  fNarip,p)  dp)- Cipf\  (^-39) 

The  rest  of  this  argument  justifying  (5.14)  is  routine  but  messy  and  so  we  transfer  it 
to  Appendix  A.  The  conclusion  of  our  argument  shows  that  whenever  the  conditions 
0<  T<-j^  and  (5.16),  which  is  n>  40  {2  +  Aa)j  [9  (pp‘‘  '^),  hold,  we  have  (5.14),  which 

is 

Prob(|M(p,p,n)-£((^/p)^•r  <p)|<  t)>1-<p  (5.40) 

Finally,  we  must  derive  the  bound  (5.17)  from  (5.13)  and  (5.14).  Once  this 


derivation  is  complete,  so  is  the  Justification  of  the  Statistical  Relative  Error  Bound.  It  is 
fairly  clear  what  must  be  done  to  derive  (5.17),  but  unfortunately  the  details  are  again 
messy,  so  we  place  them  in  Appendix  B. 
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6.  Estimating  the  Correlation  Dimension  for  the  Lorenz  and  Henon  Attractors 

A  fundamental  problem  in  interpreting  an  estimate  of  the  correlation  dimension  such 
as  that  given  in  Figure  1  for  the  Lorenz  attractor  is  determining  the  reliability  of  this 
estimate.  With  the  integral  method  for  estimating  d,  we  may  address  this  problem  using 
two  distinct  methods.  First,  by  evaluating  the  approximants  of  d  given  by  (5.12)  for 
many  values  of  p  and  p,  we  may  identify  those  values  of  p  for  which  these  approximants 
cluster  somewhat  as  the  value  of  p  varies.  Second,  we  may  use  the  statistical  techniques 
introduced  in  section  5  to  find  the  probability  that  the  numerically  estimated  value  of  d^  lies 
with  probability  at  least  l-<p  within  a  calculable  error  bound  £(Ai,  t,p,dj  (5.18)  of  the 
actual  value  d. 

We  may  illustrate  the  first  method  by  estimating  the  value  of  the  correlation 
dimension  d  for  the  standard  Lorenz  and  Henon  attractors  (with  parameter  values  Ra  =  2S, 
<T=10,  fc  =  |  as  in  Lorenz,  1963;  and  a„  =  1.40,  b„  =0.3  as  in  Henon,  1976).  For  each 
attractor,  we  use  the  values  p  =  ^,-j,|,l,|, 2,3,4  to  estimate  the  limiting  moment  M(p) 
using  the  approximate  moment  M{p,p,n)  (5.11)  from  a  sequence  of  distances 

numerically  generated  in  the  usual  way.  Here  we  use  the  classical  method  for  generating 
distances  rj,r2,r3,...  in  the  order  they  are  calculated  from  a  single  generic  trajectory;  as 
noted  in  section  5,  this  is  essentially  equivalent  to  the  method  introduced  in  that  section  for 
generating  the  distances  from  a  generic  pair  of  trajectories. 

In  Figure  3,  we  plot  the  curves  d=-d^  =d(p,p,n)  given  by  (5.12)  for  varying  values 
of  p  for  the  Lorenz  attractor  (Figure  3a;  in  which  the  number  k  of  transient-free  time  series 
points  is  20,000)  and  the  Henon  attractor  (Figure  3b;  ^  =  6000);  here  the  number  of 
distances  n  =  k{k-l)/2  is  given  by  (5.4).  As  we  note  below,  only  a  fraction  of  these 
distances  contribute  to  the  figures.  Using  any  one  of  these  curves,  we  would  look 
classically  for  a  scaling  region  in  which  the  curve  is  relatively  constant,  and  we  would 
declare  this  constant  value  to  be  our  estimate  for  the  correlation  dimension.  As  is  fairly 
clear,  such  a  flat  segment  might  be  found  between  approximately  p=1.2  and  p=1.7 
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Estimates  of  d  for  8  values  of  p 
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Figure  3.  Eight  estimates  d^=d{p,p,n)  (5.12)  of  the  correlation  dimension  d  for  the 
standard  Lorenz  (Ra=2S,  <7=10,  (a)  and  Henon  {a„  =1.40,  b„  =0.3)  (b) 

attractors  using  p  =  ^,y,|,l,|,2,3,4;  in  both  cases,  the  curves  vary  monotonically 
from  the  one  for  p  =  |  to  the  one  for  p  =  4.  For  the  Lorenz  attractor,  20,000 
transient-free  time  series  points  and  10,000  bins  are  used,  while  for  the  Henon 
attractor  6000  points  and  200  bins  are  used,  although  only  a  fraction  of  these 
contribute  to  the  given  ranges  in  p. 
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for  the  Lorenz  attractor  (Figure  3a),  but  it  is  impossible  to  find  for  the  Henon  attractor 
(Figure  3b).  The  relatively  flat  segment  for  the  Lorenz  attractor  is  also  a  region  in  which 
the  dimension  curves  cluster,  giving  values  for  d  between  2.08  and  2.09;  these  values  are 
close  to  those  given  by  the  slope  estimates  for  small  values  of  p  in  Figure  la  as  well  as  to 
those  in  the  literature  (e.g.,  Grassberger  and  Procaccia,  1983a,  b  and  Nese  et  al.,  1987 
report  d  =  2.05 ±0.01;  Simm  et  al,  1987  find  t/  =  2.16±0.3).  In  contrast,  the  estimates  of 
d  for  the  Henon  attractor  do  not  cluster  consistently  over  any  range  of  p,  leading  to  a  large 
uncertainty  in  the  estimate  for  d.  The  best  we  can  say  is  that  d  apparently  lies  between  1.17 
and  1.27  for  p  between  0  and  0.5  or  between  1.21  and  1.27  for  p  between  0.25  and  0.45; 
these  values  are  nevertheless  in  the  range  reported  in  the  literature  (e.g.,  Grassberger  and 
Procaccia,  1983a,  b  cite  1.21  ±0.01  when  the  full  map  is  used,  but  1.25 ±0.02  when  the 
time  series  from  a  single  variable  is  used;  Abraham  et  al.,  1986  find  1.24±0.02  when 
10,000  points  are  used;  Simm  et  al,  1987  report  1.25±0.1  when  8192  points  are  used; 
Ameodo  etai,  1987  obtain  1.199±0.003  when  10*  points  are  used  with  p  as  small  as  2'^'; 
and  Grassberger,  1988  finds  1.2±0.05  when  4x10^  points  are  used  with  p  as  small  as 
2"^^).  Nevertheless,  Theiler  (1988)  notes  that  the  Henon  attractor  is  probably  nonlacunar, 
which  in  our  context  implies  that  a  =  A  (e.g.,  Grassberger,  1988),  and  therefore  that  the 
limiting  moment  M{p)  exists,  so  that  our  integral  method  results  probably  can  be  trusted. 

By  way  of  caution,  we  recall  that,  as  remarked  at  the  beginning  of  section  5, 
conventional  error  estimates  lead  to  predicted  errors  so  unrealistically  small  that  at  least  one 
of  the  following  three  possibilities  must  hold  in  order  to  invalidate  these  estimates: 

i)  The  number  n  is  not  large  enough  to  justify  an  appeal  to  the  Central 
Limit  Theorem, 

ii)  The  observations  are  not  sufficiently  independent,  or 

iii)  The  values  of  r  used  are  not  sufficiently  small  for  E{{rlpY  \r<(Yj  and 
M(p)  to  be  interchangeable. 
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It  seems  to  us  that  it  is  the  third  of  these  that  is  the  likely  culprit.  Thus,  as  in  section  5,  we 
turn  from  the  conventional  approach  to  one  yielding  coarser  but  more  reliable  bounds,  using 
our  hypothesis  that  C(r)  ~  r'‘  holds  in  the  strong  sense  to  diminish  some  of  the  uncertainty 
introduced  by  possibility  iii)  above. 

To  apply  this  second  approach,  we  need  to  associate  the  probability  1  -  with  the 
number  n(p,  t,  <p)  of  independent  distances  used  to  estimate  d  in  the  empirical  scaling 
region  of  p.  As  noted  above,  this  scaling  region  is  given  classically  by  0  < p^„  < p< Po,  in 
which  p^„  is  the  smallest  distance  for  which  the  dimension  curves  are  relatively  level  and 
Po  is  the  largest.  This  is  also  a  reasonable  way  to  choose  the  scaling  region  in  our  case.  To 
choose  appropriate  values  of  Po,  we  first  note  that  as  the  value  of  p  increases,  the  function 
A#(p,  p,  n)  approaches  0  for  any  value  of  p  and  n.  Thus  it  immediately  follows  that  the 
function  pM(p,p,n)/{\-M(p,p,n))  in  (5.12)  approaches  the  value  of  0  for  large  values  of 
p.  Consequently,  we  expect  the  beginning  of  a  plunge  in  the  graph  of 
pM{p,p,n)/{l- M(,p,p,n))  to  signal  the  transition  beyond  the  scaling  region  of  p.  Thus, 
we  choose  the  largest  value  Po  of  p  to  be  the  one  that  marks  the  beginning  of  the  steady 
decrease  to  0  in  the  graphs  of  pM{p,p,n)f{l- M{p,p,n))  such  as  those  in  Figures  3a  and 
3b.  Next  we  choose  a  value  for  p^„  <  p^  that  is  large  enough  to  exclude  the  noisy, 
undersampled  region  for  small  values  of  p. 

Once  we  have  determined  the  empirical  scaling  region,  we  may  specify  the  other 
parameters  needed  to  estimate  the  relative  error  bound  for  d.  A  range  of  trial  dimensions  d^ 
that  approximate  d^  is  given  by  the  variation  within  this  scaling  region  of  d^  with  p.  As 
discussed  further  below,  the  values  of  a  and  A  are  given  by  the  minimum  and  maximum 
values  of  the  ratio  C(p)lp‘‘^  within  the  scaling  region.  By  definition,  we  have 
n(p,  <p,  T)  -  nC(Po ) .  Thus,  after  reexpressing  (5. 16)  as 

,  40(2/4“ +(2/4) 

— - 

9a-p^nC{p,) 


(6.1) 
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we  are  able  to  find  acceptable  values  for  <p  and  t.  If  Po  <  1,  then  we  choose  the  value  of  dj 
in  (6.1)  to  be  the  upper  bound  of  the  range  for  d,-,  while  if  Po  >  1,  then  we  choose  the  lower 
bound.  We  may  now  use  (5.18)  to  obtain  an  initial  estimate  for  the  relative  error  bound 
eiAi,r, p,dj  for  d  with  the  associated  probability  l-<p  and  moment  error  tolerance  t. 
This  value,  however,  is  not  the  smallest  one  available  because  e(Aj,  T,p,d^)  depends  on 
both  p  and  d^,  which  are  parameters  that  may  be  regarded  as  varying  independently  within 
suitable  bounds. 

To  find  this  minimum  value  of  the  relative  error  bound,  we  thus  regard  as  fixed  the 
values  of  rand  Ap  which  is  given  in  (5.19)  by  the  known  values  of  a  and  A.  Then  we  hold 
d^  constant  and  minimize  the  value  of  e(A,,  r,p,dj  with  respect  to  p.  Using  the  symbolic 
manipulator  Derive,  we  find  that  this  minimum  value  is  given  by 

[V(1-A,)->/t(A,  +  t)][Vt(A,  +  t)(1-A,)  +  A,-i] 


We  note  that  this  minimum  value  is  actually  independent  of  d^  and,  because  it  is  attained 
for  every  value  of  d^,  it  must  be  the  global  minimum  Cn^JA,,  T)  =  e^„(A,,  x,d^)  that  we 
seek.  The  value  of  p  yielding  (6.2)  is  given  by 


P  min 


_ 


T+1 


fT(l-A,) 
A,  +  t 


-  T 


for  A,  <  1 


(6.3) 


There  are  two  limiting  cases  of  interest.  First,  for  the  case  t  =  0  of  perfect  estimates 
for  E({rlpY:r<p),  we  cannot  find  a  relative  error  bound  smaller  than 


e_(A,,0)  = 


l-A, 


(6.4) 


for  any  t.  We  call  this  the  geometric  part  of  the  error  because  it  is  intrinsic  to  the  geometr>' 
of  C(p),  unlike  the  error  t  that  is  due  to  noise  introduced  by  undersampling  at  small  values 
of  p.  In  this  perfect  case,  must  vanish,  which  is  forbidden  because  p  >  0  is  required  in 
(5.18).  However,  we  may  interpret  =  0  as  meaning  that,  no  matter  what  value  of  p  we 
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use  to  estimate  a  relative  error  bound  for  d,  we  could  find  a  better  estimate  using  a  still 
smaller  value  of  p.  Second,  for  the  case  A,  =  0,  we  obtain  the  limiting  relative  error  bound 

eau.(0.-r)  =  77^  (6.5) 

(l-T) 

with  corresponding  given  by 

(6.6) 

1+  T 

We  call  e„ua(0,  t)  the  statistical  part  of  the  error  because  it  arises  from  random  errors  in 
measurement.  Below  we  calculate  (6.2),  (6.4)  and  (6.5)  to  see  how  close  our  error  bound 
Cmin(^]>  is  to  the  limiting  error  bounds  e„j„(Ai,  0)  and  e^„(0,  t). 

For  the  Lorenz  attractor,  we  inspect  Figure  3a  to  obtain  the  estimates  Po  -1.'/  and 
Pmin  for  the  scaling  region,  while  for  the  Henon  attractor  we  obtain  -0.45  and 
Pmin  -0.25.  The  value  of  nip,  <p,  t)  ~nC(Po)  for  the  Lorenz  attractor  is  nip,  (p,  t)  - 1  x  10^ 
while  for  the  Henon  attractor  it  is  nip, <p,t)- 5 x  10^ .  To  find  the  values  of  a  and  A,  we 
calculate  the  ratio  C(p)/ p'^^  for  several  trial  dimensions  that  are  suggested  by  the  plots 
of  d^  versus  p  in  Figure  3.  It  is  clear  from  these  figures  that  2.08  <  <4  ^  2.10  for  the  Lorenz 
attractor,  while  1.21  <  <4  ^  1-27  for  the  Henon  atu-actor.  Within  the  empirical  scaling  region 
Pmin  <  P<  Po>  ihe  indicated  value  for  a  is  the  minimum  value  of  the  ratio  Cip)lp‘‘^ ,  while  A 
is  the  maximum  value,  as  given  by  (2.7).  Of  course,  by  definition  a  and  A  are  the  minimum 
and  maximum  values  of  Cip)lp‘‘^  over  the  entire  range  0<p<Po.  But,  owing  to  the 
unavoidable  undersampling  of  the  attractor  given  by  our  finite  dataset,  we  must  attempt  to 
estimate  the  values  of  a  and  A  from  the  behavior  of  Cip)fp‘^^  within  only  the  chosen 
scaling  region.  Finally,  we  note  that,  if  Theiler  (1988)  is  correct,  then  the  actual  values  of  a 
and  A  should  be  equal  for  the  Henon  attractor  (Grassberger,  1988).  However,  we  will  not 
assume  them  to  be  equal  in  our  error  estimates  below,  and  so  use  the  following  numerical 
values  for  a  and  A  giving  A,  ^0.  For  the  Lorenz  attractor,  we  find  a -0.00215  and 
A -0.00217  (not  shown),  and  for  the  Henon  attractor,  a -0.716  and  A -0.778  (Figure  4). 
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Figure  4.  Eslimales  of  the  ratio  C(p)lp‘‘'  for  1 1  trial  values  dj  of  the  correlation  dimension 
d  for  the  Henon  attractor  {a„  =  1.40,  b„  =0.3).  The  maximum  value  of  this  ratio  in 
the  empirical  scaling  region  p^,„<p<P(,,  given  by  Pm,„=0.25  and  o„=().45  and 
denoted  by  the  dashed  lines,  provides  an  estimate  for  A,  while  the  minimum  value 
gives  a.  Here  6000  points  and  200  bins  are  used,  although  only  a  fraction  of  the.se 
contribute  to  the  distance  interval  used. 
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Our  next  task  is  to  estimate  the  product  (pT^  using  (6. 1 ),  under  the  assumption  that 
we  actually  have  nip,  r,  <p)  independent  evaluations  of  the  distances.  With  the  above  values 
of  a.  A,  dj,  and  nip,  x,(p),  we  find  that  <pT^  -2x10"^  for  the  Lorenz  attractor,  while 
<px^  -  lx  10”’  for  the  Henon  attractor.  Estimates  using  (6.2)  show  that  acceptable  relative 
error  bounds  of  no  more  than  15%  to  20%  are  only  achieved  if  t<0.03.  With  <ip^0.2  to 
give  a  probability  of  at  least  80%  that  we  have  a  good  estimate  for  d,  we  conclude  that  we 
would  need  to  have  ^  2x  10^.  Although  we  used  20,000  points  giving  approximately 
2x10*  distances  in  our  estimates  for  the  Lorenz  attractor,  we  have  only  10%  of  the 
distances  we  need  in  the  scaling  region  <  p<Po  for  us  to  find  acceptable  relative  error 
bounds  for  d;  such  a  calculation  would  require  200,000  time  series  points,  which  is  10  times 
that  normally  considered  in  an  estimate  of  d  (e.g.,  Grassberger  and  Procaccia,  1983a,  b).  In 
contrast,  the  6000  points  we  used  for  the  Henon  attractor  give  enough  distances  in  the 
scaling  region  to  allow  us  to  produce  acceptable  relative  error  bounds  for  d. 

We  now  are  ready  to  estimate  the  relative  error  bounds  e(A,,T,p,J^)  (5.18), 

(^.4),  and  e^„{0,  T)  (6.5)  for  the  Henon  attractor.  Two 
combinations  of  (p  and  t satisfying  the  above  condition  on  tp'^  ~  1  x  lO”’  are  used:  t=  0.02, 
<p  =  0.025  and  t=0.01,  (p  =  0.1.  These  results  are  reported  in  Table  1  and  illustrated  in 
Figures  4  and  5.  Inspection  of  Table  1  reveals  that  the  smallest  error  occurs  for  dj  =  1.26, 
and  inspection  of  Figure  4  shows  that  the  ratio  C(p)/p‘^'^  varies  the  least  in  the  empirical 
scaling  region  for  this  value  of  dj.  The  geometric  part  of  the  error  implies  an 

estimate  of  at  least  dj  =1.25±0.02,  while  the  statistical  part  implies  d^  =  1.2510.1  with 
probability  97.5%  or  =  1-2510.05  with  probability  90%.  In  contrast,  the  complete 
relative  error  bound  e„u„(Ai.T)  leads  to  the  estimates  dj  =  \.25±0.\5  with  probability 
97.5%  and  dj  =  1.2510.1  with  probability  90%.  Incidentally,  the  second  complete  estimate 
agrees  with  the  more  ad  hoc  one  of  Simm  et  al.  (1987)  and  the  geometric  part  agrees  with 
the  ad  hoc  one  of  Grassberger  and  Procaccia  (1983a,  b). 


Table  1.  Estimates  for  the  Henon  attractor,  expressed  as  a  percent  of  d,  of  the  relative  errors  e(A,,  t,  p,d„)  when  p  =  1  (5.18), 
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P 


Figure  5,  Illustration  of  how  the  optimal  value  dj  for  the  correlation  dimension  d  can  be 
found  by  minimizing  the  variation  of  C{p)lp‘‘^  for  the  Henon  attractor  {a„  =1.40, 
bf,  =  0.3).  The  curve  dj  =  \.2(>  corresponds  to  the  case  in  Table  1  that  minimizes  the 
height  A- a  of  the  shaded  rectangle  in  the  figure.  The  empirical  scaling  region 
Pmin  <  P  <  Po  is  indicated  by  the  dotted  lines. 


93 


Although  our  estimates  are  very  coarse,  they  are  objectively  determined  and  identify 
the  most  likely  value  for  the  correlation  dimension  d  that  can  be  obtained  from  a  transient- 
free  time  series.  Inspection  of  (6.2)  and  Table  1  reveals  that  the  values  for  are  not 
needed  to  identify  the  optimum  value  of  d.  The  error  bounds  are  closest  when  the  value  of 
A;  is  smallest,  and  from  its  definition  (5.19),  we  see  that  this  occurs  when  the  difference 
A -a  is  smallest.  Thus,  we  need  only  plot  the  ratio  C{p)lp‘‘^  as  a  function  of  p  for  various 
candidate  values  dj,  and  then  seek  the  value  producing  the  smallest  difference  A -a. 
These  candidate  values  are  given  by  the  approximants  d^  that  are  calculated  using  the 
moment  M{p,p,n)  via  (5.12). 

The  procedure  of  minimizing  A-a  h  illustrated  in  Figure  5,  in  which  we  reproduce 
one  of  the  curves  in  Figure  4.  As  indicated  by  the  more  thorough  analysis  in  Table  1,  the 
difference  A-a  is  smallest  for  the  case  dj  =  \.26  for  the  Henon  attractor.  Although 
finding  the  values  of  a  and  A  that  minimize  A-a  leads  to  the  geometric  part  e„^„(Ai,())  of 
the  relative  error  bound  for  d,  it  does  not  reveal  anything  about  the  statistical,  and  major, 
part  resulting  from  the  moment  error  tolerance  t  (see  Table  1).  Nor  does  it  reveal  the 
magnitude  e(A,,  T,p,dJ  or  e^„(Ai,  t)  of  this  error.  We  therefore  recommend  that  a  full 
analysis  of  the  relative  error  bounds  be  performed  to  see  if  they  are  small  enough  to  allow 
acceptance,  with  reasonably  high  probability,  of  the  estimated  values  for  d. 


7.  Hentschel-Procaccia  Dimensions 

Hentschel  and  Procaccia  (1983)  introduce  a  scale  of  dimensions  generalizing  the 
Grassberger-Procaccia  (1983a,  b)  correlation  dimension  d.  The  same  device  that  we  have 
used  in  section  2  may  be  used  to  improve  the  numerical  estimates  of  the  Hentschel  and 
Procaccia  dimensions  as  well.  However,  we  first  need  to  establish  some  background  and 
notation  before  applying  tha'  device. 
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We  assume  that  we  have  a  C‘  diffecmorphism  F:R^ and  a  probability 
measure  /i  such  that  Fx  •••  xF:  x  x  is  ergodic  with  respect  to 

/tx-"X/r.  Then  for  <3  =  1,2, 3,...  we  may  define  the  corresponding  HenLschel-Procaccia 
dimension  D(q)  by  setting  first 


C^(r)  =  j/i(5U,r))’tf/iU)  (7.1) 

where  B(x,  r)  is  the  ball  of  radius  r  and  center  x, 

fi(jc,r)  =  {yl|jc-y|<;}  (7.2) 

Then  we  define  with  the  limit  equation 

1  In  C.(r) 

D,=-lim— ^  (7.3) 

q  '-’0  Inr 

provided  this  limit  exists.  We  recall  from  (2.2)  that  the  correlation  dimension  d  =  D^. 

We  note  that  we  have 


^{B{x,r))  =  jx(\x-y\/r)d^i(y)  (7.4) 

where  x(^)  is  given  by  (5.10),  so  that 

C^ir)  =  l  ■  ■■\ Xi\x-yi  \/r)  ■  ■  ■  x{\x-}\ \/r)dm}\  )■■■  d^{y^)d^2{x)  (7.5) 

holds.  We  define  a  function  r^{x,y,,...,y^)  by  setting 

r,(x,y,,....y,)  =  min{ljf-yi!,...,|x->J}  (7.6) 

Then  we  have 

C,(r)  =  J  •  ■  •  j x[r^(x, y, , . ..,y^)lr)d^{\\  )  ■■■  d^i{y^)d^{x)  (7.7) 

Now  we  introduce  a  family  of  limiting  moments  M^{p)  by  .setting 

Af/p)  =  lim£((rjp)'’:r,<p)  (7.8) 


provided  that  the  limit  exists.  We  note  that  we  have  Mj  (p)  =  M(p). 
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where  M,  {p,  y)  =  M{p,  y). 
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8.  Concluding  Remarks 

In  this  article  we  present  a  new  method  for  estimating  the  value  of  the  Grassberger 
and  Procaccia  (1983a,  b)  correlation  dimension  d.  This  method  connects  the  dimension  at 
the  origin  of  a  continuous  probability  distribution  C(r)  on  [0,  <»)  with  its  limiting  moment 
M{p)  at  the  origin.  More  specifically,  we  define  the  dimension  of  C(r)  at  the  origin  to 
have  the  value  of  d  if  the  relation  C(r)  ~  r‘‘  holds  in  the  weak  sense  (defined  in  section  2)  as 
r  approaches  zero.  If  the  dimension  has  a  finite  value  d,  then  the  expected  value  of  the 
random  variable  {rlpY,  relative  to  the  condition  r  < p,  has  limit  M(p),  given  by  (2.5),  as  r 
approaches  zero.  The  connection  between  the  limiting  moment  and  the  dimension  d  is 
given  by  the  simple  formula  (2.6), 


pMjp) 

l-M(p) 


(8.1) 


The  classical  algorithm  for  estimating  the  dimension  identifies  it  with  the  limit  given 
by 


'^-♦0  dlnr  '-*°C{r)  dr 


(8.2) 


Therefore,  this  algorithm  requires  the  function  C(r)  to  be  differentiable,  and  then  the  slopes 
of  the  function  y  =  InC(expjc)  approach  d  as  x  approaches  -oo  (where  x  =  Inr).  However, 
because  numerical  differentiation  is  very  sensitive  to  errors  (see  Figure  la),  various  kinds  of 
mean  slopes  are  used  instead  to  generate  a  numerically  stable  algorithm. 

The  method  we  introduce  here,  leading  to  (8.1),  uses  the  hypothesis  that  C(r)  is 
continuous.  Instead  of  differentiating,  we  integrate  to  define  the  function  /(p)  given  by 

I{p)  =  j'C{r)dr  (8.3) 

We  easily  see  that  this  new  function  satisfies  the  relation  /(p)  -p*'*'  so  that 
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P-*o  /(p) 


(8.4) 


holds,  provided  the  limit  exists.  This  formula  is  similar  to  one  suggested  by  Theiler  (1988, 
1990). 

By  using  (8.4)  in  place  of  (8.2),  we  eliminate  the  need  for  numerical  slope-finding 
routines,  replacing  them  with  the  more  stable  numerical  integration  routines  required  for 
tabulating  the  new  function  /(p)  (see,  for  example.  Figure  lb).  In  addition,  by  using 
algorithms  based  on  (8.4)  instead  of  (8.2),  we  also  extend  the  applicability  of  the 
dimension-finding  algorithms  to  the  important  case  for  which  C(r)  is  not  differentiable. 
However,  to  improve  directly  the  accuracy  of  an  algorithm  based  on  (8.4),  we  require  more 
accurate  numerical  integration  of  the  nondifferentiable  monotone  continuous  function  C(r). 
Because  an  algorithm  other  than  Riemann  sums  is  not  available,  we  use  instead  a  Stieltjes 
integration  by  parts  to  transform  (8.4)  into  (8.1),  which  requires  only  evaluation  of  certain 
means,  thereby  bypassing  errors  introduced  by  using  inadequate  numerical  integration 
routines.  By  evaluating  these  means,  we  eliminate  the  need  to  order  by  magnitude  the  set  of 
distances,  as  is  often  done  in  the  classical  slope  procedure  (Albano  et  ai,  1988)  but  that  is 
computationally  very  expensive  for  datasets  of  the  large  size  required  to  produce  usable 
error  bounds  (see  section  6).  Finally,  integral  methods,  based  on  (8.2)  or  (8.1),  make  it 
unnecessary  to  find  the  function  C(r)  before  finding  the  dimension  d. 

A  further  advantage  of  integral  methods  over  slope  methods  is  the  apparently 
superior  rate  of  convergence  of  the  former  (Theiler,  1988).  In  any  numerical  or  empirical 
approximation  of  the  probability  distribution  C(r),  the  smaller  distances  are  noise- 
dominated  in  comparison  with  the  larger.  What  results  is  a  function  C,(r)  approximating 
C(r)  well  for  large  values  of  r  and  poorly  for  small  values  of  r.  As  Ben-Mizrachi  et  al. 
(1984),  Simm  et  al.  (1987),  and  Theiler  (1986,  1991)  note,  C,(r)  scales  with  the  dimension 
d  for  large  values  of  r,  but  with  the  much  larger  embedding  dimension  c  -  2d-t- 1  for  small 
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values  of  r  when  the  correlation  integral  estimate  is  contaminated  by  noise.  We  model  this 
situation  by  considering  the  special  nonlacunar  case  of  an  actual  distribution  C(r)  =  that 
is  approximated  by  a  continuous  distribution  C,(r)  given  by 


C,{r)  = 


for  0<r<5 
for  5<r<p 


(8.5) 


in  which  5  is  a  positive  number.  Thus  we  regard  those  values  of  r  lying  below  5  as 
belonging  to  the  noise-dominated  region  and  those  lying  above  5  as  belonging  to  the  well- 
sampled  region,  with  the  transition  from  one  region  to  the  other  taking  place  abruptly  at 
r  =  5.  Then  we  compare  the  slope  method  with  the  integral  method  by  using  the  same 
working  interval  [5,,p],  with  0<  5,  <  5<p<  1,  to  fit  a  line  in  the  lnC(r)-lnr  plane  for 
the  slope  method  and  to  integrate  C^(r)  for  the  integral  method.  Thus  the  first  procedure 
produces  the  mean  slope  over  the  interval  [5,,p]  and  the  second  produces  an  approximation 
of  /(p).  From  each  of  these  results,  we  derive  the  corresponding  estimate  of  d.  Upon 
comparing  these  two  estimates  as  the  proportion  {5-5^)lp  of  the  noise-dominated  region 
decreases  to  zero  through  positive  values,  we  find  analytically  in  section  3  that  the  integral 
estimate  approaches  the  dimension  d  much  more  rapidly  than  does  the  slope  estimate.  In 
addition,  we  find  numerically  an  even  more  dramatic  improvement  in  the  rate  of 
convergence  (see  Figures  2b,  c),  strongly  suggesting  that  the  integral  method  is  far  less 
sensitive  to  contamination  by  errors  at  small  values  of  r  than  is  the  slope  method.  Of 
course,  considerable  work  remains  to  be  done  in  this  direction:  first,  to  compare  these  two 
methods  in  paradigmatic  cases  analogous  to  those  given  by  (8.5),  but  with  a  gradual 
transition  from  noise-contaminated  to  well-sampled  regions,  and,  second,  to  compare  the 
two  in  the  case  of  a  similar  blend  of  more  general,  especially  lacunar,  distributions  (e.g., 
Theiler,  1988). 
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Furthermore,  in  section  7  we  extend  the  integral  method  to  apply  to  the  case  of  the 
HenLschel-Procaccia  generalized  dimensions  of  integral  index  (Hentschel  and  Procaccia, 
1983),  and  in  section  4  we  also  recover  the  formulas  of  Takens  (1985). 

Besides  the  apparently  more  rapid  convergence  of  the  integral  method  over  the  slope 
method,  a  second  advantages  lies  in  its  formulation  in  terms  of  expectations  and 
probabilities;  This  formulation  enables  us  to  use  the  simple  Kolmogorov  Inequality  to 
bound,  in  an  a  priori  manner,  the  probability  that  our  estimate  of  the  dimension  differs  from 
the  actual  dimension  by  more  than  a  prescribed  amount.  In  this  article  we  begin  this 
process  by  arriving  at  somewhat  coarse  error  bounds,  and  identifying  four  obstacles  to 
obtaining  more  refined  and  satisfactory  versions  of  these.  First,  our  results  depend  upon  the 
hypothesis  that  the  relation  C(r)~r‘‘  holds  in  the  strong  sense  that  is  defined  in  section  2; 
however,  it  is  fairly  clear  how  to  adjust  our  arguments  to  the  weak  case  so  that  this  obstacle 
is  not  a  very  serious  one.  Second,  our  arguments  require  that  two  successive  evaluations  of 
certain  random  variables  be  nearly  independent.  This  requirement  may  be  met  by  thinning 
out  the  time  series,  by,  for  example,  using  every  tenth  point  encountered  to  calculate  the 
series  of  distances  used  to  estimate  C(r).  Unfortunately,  such  a  thinning,  together  with  the 
required  length  of  the  surviving  time  series,  amounts  to  restricting  the  validity  of  the  error 
bounds  we  find  to  cases  having  extremely  long  time  series,  such  as  computer-modeled 
attractors.  The  third  obstacle  results  from  the  need  to  bound  the  magnitude  of 
^E{(r/py  :r<  p^- M(p)  .  In  this  article  we  make  use  of  the  very  weak  bound 
A,  =  Aja-ajA,  which,  in  addition,  makes  sense  only  when  the  relation  C{r)~r‘‘  holds  in 
the  strong  sense  for  which  a  ^  C{r)lr‘‘  <  A .  What  we  need  here  is  a  stronger  bound, 
especially  one  that  is  valid  in  the  weak  case.  Finally,  the  fourth  and  most  serious  obstacle  is 
posed  by  those  correlation  integrals  C{r)  for  which  the  limiting  moment  M{p)  does  not 
exist  (Theiler,  1988).  As  noted  in  section  2,  the  limiting  moment  along  an  appropriate 
subsequence  p,,p2,pj,...  still  exists  and  yields  the  correlation  dimension  via  (2.21),  but  the 
problem  of  extracting  this  subsequence  is  formidable. 
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An  alternative  approach  would  be  to  develop  usable  algorithms  to  enable  us  to 
decide  when  a  specific  attractor  yields  a  strong  case  of  the  relation  C(r)~r‘',  not  only  to 
detect  the  strong  case,  but  also  to  provide  estimates  of  the  associated  parameters  p^,  A  and 
a.  For  example,  we  may  ask  whether,  for  all  practical  cases,  a  numerical  model  of  the 
Lorenz  system  (Lorenz,  1963)  yields  a  strong  case  of  the  relation  C(r)  -  r‘‘ .  Or,  even  more 
specifically,  whether  we  may  prove  at  least  that  the  Smale-Williams  hyperbolic  attractor 
(Smale,  1977)  yields  a  strong  case,  and  if  so,  whether  we  may  estimate  rigorously  values  for 
the  parameters  Po ,  A  and  a.  We  suggest  that  such  considerations  as  these  may  lead  to  an 
objective  and  effective  algorithm  to  determine  the  scaling  region  in  a  computational  or 
observational  setting. 

Application  of  our  coarse  relative  error  bound  algorithm  to  the  standard  Henon 
attractor  (Henon,  1976)  yields  values  of  the  correlation  dimension  equal  to  dj  =  1.25±'J.15 
with  probability  97.5%  and  moment  error  tolerance  t=0.02,  and  d;^  =  1.25±0.1  with 
probability  90%  and  t=0.01,  respectively.  These  error  bounds  are  similar  to  those 
reported  in  the  literature  (e.g.,  Simm  et  al,  1987).  Use  of  longer  series  to  reduce  noise 
contamination  at  small  distances  would  presumably  lead  to  improved  estimates.  If  we  are 
not  interested  in  error  bounds  on  the  value  of  the  dimension  d,  then  our  analysis  suggests  a 
simple  procedure  for  estimating  d  from  a  time  series  of  measurements.  First  we  find  several 
independent  estimates  d^  for  d  using  the  limiting  moment  M{p).  Upon  plotting  the 
approximate  values  of  as  functions  of  the  maximum  distance  p,  we  obtain  a  relatively 
small  range  of  candidates  dj  for  in  an  appropriate  working  interval  of  p.  We  determine 
the  optimal  value  of  d  from  this  set  of  candidates  by  determining  which  one  yields  the 
smallest  difference  A-a  in  the  ratio  C{p)lp'‘^  over  the  associated  working  interval  of  p. 
This  approach  is  used  by  Fosmire  (1993)  in  an  analysis  of  two  time  series  of  atmospheric 
boundary  layer  winds. 
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The  probabilistic  context  associated  with  our  derivation  leads  to  another  method  for 
generating  the  correlation  function  C(r).  Instead  of  using  all  the  mutual  distances  between 
points  on  a  relatively  short  time  series  of  points  on  the  attractor,  we  may  use  the  distances 
between  corresponding  points  on  two  extremely  long  time  series.  Because  the  two  long 
time  series  pass  through  more  of  the  attractor  at  a  greater  variety  of  scales,  we  may  expect 
that  the  approximation  to  C(r)  thereby  obtained  is  more  accurate  than  the  one  obtained 
from  the  mutual  distances  in  a  short  series.  Of  course  the  principal  drawback  to  this 
technique  is  that  it  requires  such  a  long  time  series  that  the  method  is  effectively  restricted 
to  computer-modeled  attractors. 
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Appendix  A 

Completion  of  the  Proof  of  the  Statistical  Error  Estimate  (5.14) 

Our  next  task  is  to  find  an  upper  bound  for  the  values  of  Var(p,  p)  and 
C(p)-C(pf ;  we  may  combine  (5.23)  and  (5.29)  with  the  definition  (5.32)  to  obtain 


Var(p,p)< 


A  2p  a  2p  \ 

<’  Tl 

a  p+d  A  2p+d 

^p+d) 

C(p) 


(Al) 


We  note  in  passing  that  the  left  side  of  (Al)  is  indeed  positive.  Using  (Al),  we  may 
weaken  (5.39)  to  the  inequality 


Prob(|A/o(p,p,/i)-£((r/p)'’;t(r/p))|<p‘'^*  and  |  A(p,n)/n-C(p)|  <  p'^^*) 


C(p) 

A  2p 

A 

2t) 

a  p+d 

p‘‘*^^n 

p2d*2S 

(A2) 


We  check  that  the  two  premises 

Mo  (p,  p,  n)  -  E[{rlpy  x{rlp))\  <  P' 


d*S 


and 


\N(p,n)/n-Cip)\<  p' 


d*S 


(A3) 


(A4) 


imply  the  consecutive  inequalities 


Mo(p,p,n)  E[{rlpy  x{r!p)]  ^  p^*^  ^  E[{rjpy  xirlP))p'*^ 

Nip,n)/n  C(p)  N{p,n)!n  C{p)N{p,n)/n 

^  ^  Ejjrlpy  x[rlp))p‘‘*^ 

~C{p)-p‘‘*^'^  {C{p)-p‘‘*^)C{p) 


(A5) 


where  we  have  used  (A4)  in  the  second  step.  We  note  that  the  two  equations 
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Mip,p,n)  = 


MQ{p,p,n) 

N(p,p)/n 


r.,,  ,  ,  1  E(MpYx(rlp)) 

E{(rlpY.r^ph - - 

follow  from  their  definitions;  further,  we  note  that  the  trivial  inequalities 

^({rlpYxirlp))  ^  C(p) 
C{p)-p‘‘*^>{a~p^)p‘‘ 

hold.  Then  we  may  use  (A5)-(A7)  to  conclude  that  we  have 


\M{p,p,n)-E({r/pY:r<p)\<2-^—j  (A8) 

'  '  a-p 

Thus  the  conjunction  of  (A3)  and  (A4)  implies  (A8),  and  so  the  event  described  by  (A3) 
and  (A4)  is  contained  in  the  event  described  by  (A8).  Consequently,  we  have 


Prob|^|  M(p,  p,  n)-E[{rlpY :  r  <  p)|  <  2^^  j  ^ 
Prob(|Mo(p,p,r?)-£((r/p)'’;f(r/p))|<p‘'^^  and  l/V(p,/7)//7-C(p)|  <p''^^) 
By  combining  (A9)  and  (A2),  we  obtain 

Prob(^lM(p,p,.)-£((r/p)'’:rSp)|<2^]4I-^(^^^+A]-^ 


(AlO) 


For  values  of  p*  satisfying  0  <p*  <i^21,  we  have 


2  P'  ^2.1  , 

a-p*  a 


(All) 
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Making  an  obvious  simplification  on  the  right  side  of  (AlO),  we  then  obtain,  for  small 
values  of  p, 

Prob  (^1  A/(p.  p,  n)  -  £((r/p)' :  r  S  p)|  S  j  >  1  -  ^  (A  1 2) 

Finally,  we  may  rewrite  (A  12)  using  the  abbreviation  't=l.\p^  j a  to  obtain 

Prob  (I  M ip,  p,  n)  -  E({r/pY :  r  <  p)|  <  t)  >  1  -  y  ( A 1 3 ) 

that  is  valid  for  0  ^  t< Equivalently,  we  have 

Prob  (I  M (p,  p,  n)  -  E[{rlpY :  r  $  p)|  >  t)  <  y  +  -^5iE  ( ^  1 4) 

Thus,  we  arrive  at  the  following  probability  estimate: 

Statistical  Error  Estimate  If  the  conditions  0<t^^  and  (5.16),  which  is 
/i>  40(2A^  +  Afl)/(9<pp‘^aV),  hold,  then  we  have 

Prob{|  M(p,  p,  n)  -  E({r/pY :  r  <  p)|  <  t)  >  1  -  9  (A1 5) 

or,  more  precisely 

Prob(|  M  ip,  p,  n)  -  E({r/pY  :r<p)|<T)>l-9-y  ^57^  ^6) 

where  rj  is  the  deviation  from  independence  in  (5.31). 


Appendix  B 

Completion  of  the  Proof  of  the  Dimension  Relative  Error  Bound  (5.17) 


We  begin  by  noting  that,  from  (2.18)  we  have  M(p)  =  df{p+d),  so  that  we  may 
rewrite  (5.13)  as 


E[{r/pY  :r<p)- M{p)  <  A,  (l -  M{p)) 


(Bl) 


where  we  have  introduced  the  abbreviation  (5.19).  Using  (Bl),  we  see,  first,  that 


M{p)> 


^£((r/pf^r£p)-^ 
1-A, 


(B2) 


and  then  that 


I  £((r/p)'’ :  r  £  p)  -  M{p)\  ^  "  E({rlpy :  r  <  p)]  (B3) 

hold.  From  (A16),  we  have  that 

|A/(p,p,n)-£:((r/p)'’;r  <p)|<  T  (B4) 

holds  with  probability  greater  than  or  equal  to  \-^-%Qr\l9p^Wa^  provided  that  (5.16) 
holds.  Then  we  obtain  from  (B3)  and  (B4), 

\M{p,p,n)-M{p)\<-^^{x+\-M{p,p,n))+T  (B5) 

1-A, 

which  holds  with  probability  greater  than  or  equal  to  \-^-%0r]l9p^Wa^  provided  that 
n>n{p,r,(p)  in  (5.16). 

We  next  compare  the  actual  dimension  d  given  by  (2.6)  with  the  approximate 
dimension  given  by  (5.12)  that  is  obtained  from  our  approximate  moment  M{p,p,n) 
(5. 11).  In  our  calculation  below,  we  also  use 
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\-M(p)  ^  1 

l-£((r/pr:r<p)"l-A, 


(B6) 


that  follows  immediately  from  (Bl).  Then  we  combine  (Bl),  (B5),  and  (B6)  to  obtain 


d-d„ 


A,  (1-M{p,p,n))+T 


(B7) 


(l-M(p,p,n))(A/(p,p,n)-Ai  -t) 

Next  we  solve  (5.12)  for  M(p,p,n)  and  insert  the  result  in  (B7)  to  obtain  finally  our  basic 
inequality 


d-d 

d 


^  AiP+T(p-h^J  p+d^ 

~  d,-ip+dJiA^  +  r)  p 


=  e(A,,T,p,f/J 


(B8) 


bounding  the  error  \  d-d^  \  as  a  proportion  of  the  actual  dimension  d,  where  the  expression 
on  the  right  side  is  the  bound  £iA^,z,p,d^),  (5.18),  we  seek.  We  note  that  similar 
maneuvers  lead  to  a  bound  on  the  error  as  a  proportion  of  the  known  quantity  d^. 


d-d,  ^A,p+xip+dJ  p+d, 
d,  p-T(p+rfJ  d. 


The  bound  (B8)  is  estimated  in  section  6  for  the  Henon  attractor. 


(B9) 
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